A: singular gradient Matrix at initial Parameter estimates

Question

When using the following command to estimate model parameters (bell_model):

h <- c(43.34,   35.84,  33.45,  30.94,  27.35,  21.75,  13.75,  57.37,  
48.36,  44.62,  41.05,  36.49,  29.92,  21.07,  66.65,  56.65,  52.03,  
47.75,  42.54,  35.32,  25.92,  75.56,  64.60,  59.13,  54.17,  48.35,  
40.51,  30.57)
TR <-   c(2, 2, 2, 2, 2, 2, 2, 5, 5, 5, 5, 5, 5, 5, 10, 10, 10, 10, 10, 10, 
10, 20, 20, 20, 20, 20, 20, 20)
t   <- c(120,   60, 50, 40, 30, 20, 10, 120, 60, 50, 40, 30, 20, 10, 120, 
60, 50, 40, 30, 20, 10, 120, 60, 50, 40, 30, 20, 10)
dados <- data.frame(h,TR,t)

param <- list(a1 = 0.7, a2 = 0.38, a3 = 0.38, b = 0.31, a4 = 0.39)
bell_model <- nls(h ~ ((a1*log(TR)+a2)*(a3*(t^b)-a4)*41.59), dados, start = 
param)

An error occurs:

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

What would be this mistake and how to solve this problem? From now on, thank you!

Answer 1

The mistake Error in nlsModel(formula, mf, start, wts) : singular gradient matrix at initial parameter estimates means that the gradient of the search for the best estimates for your equation is unique. that is, its determinant is equal to zero and therefore the gradient matrix is not invertible.

There are many causes for this. It can range from a numerical instability (because the representation of real numbers on computers has a limited number of decimal places) to wrong initial kicks or poor specified functions. I chose to investigate the ultimate reason to solve your problem: poor specified functions.

Are you sure the model to be adjusted is h ~ ((a1*log(TR)+a2)*(a3*(t^b)-a4)*41.59)? My experience says that this model is very complicated. It is very rare to find parameter products the way you put them. So I decided to simplify this function to see what I thought.

I did some math here simplifying your formula and arrived at the model h ~ (a1*a3*log(TR)*t^b - a1*a4*log(TR) + a2*a3*t^b - a2*a4)*41.59. So I defined

a1 = a1*a3
a2 = a1*a4
a3 = a2*a3
a4 = a2*a4

Because there is no product between parameters, it seems to me much more natural to adjust a model like this than the original model. It is more natural for me and easier for the computer, because it decreases the chance of numerical instabilities.

In the end it was like this:

param <- list(a1 = 0.7, a2 = 0.38, a3 = 0.38, b = 0.31, a4 = 0.39)
bell_model2 <- nls(h ~ (a1*log(TR)*t^b - a2*log(TR) + a3*t^b + a4)*41.59, 
dados, start = param)
summary(bell_model2)
Formula: h ~ (a1 * log(TR) * t^b - a2 * log(TR) + a3 * t^b + a4) * 41.59

Parameters:
   Estimate Std. Error t value Pr(>|t|)
a1  1.26407    1.18903   1.063    0.299
a2  1.24840    1.20588   1.035    0.311
a3  4.43015    4.15081   1.067    0.297
b   0.04676    0.03749   1.247    0.225
a4 -4.69874    4.20360  -1.118    0.275

Residual standard error: 1.003 on 23 degrees of freedom

Number of iterations to convergence: 14 
Achieved convergence tolerance: 5.089e-06

Note that I used the same initial kicks of your model and my attempt found an answer. In addition, it took 14 steps to converge, which is a fairly small value.

Just be careful: a1 for me, in my parameterization, it has a different value than a1 in your parameterization. If there is a physical meaning to this parameter, it is good to see how to do the correct mathematical transformation to get its value in the context of your problem. The same goes for a2, a3 and a4.