R: singular gradient matrix at initial parameter estimates

1

When using the following command to estimate the 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 has occurred:

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

What would be this error and how to solve this problem? Thanks in advance!

    
asked by anonymous 12.05.2017 / 14:47

1 answer

0

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

The causes for this are several. It can range from numeric instability (since the representation of real numbers on computers has a limited number of decimal places) to erroneous initial kicks or poorly specified functions. I've chosen to investigate the last reason for solving your problem: poorly specified functions.

Are you sure that the template 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 it. So I decided to simplify this function to see what I thought.

I made a few accounts here to simplify your formula and I came up with the h ~ (a1*a3*log(TR)*t^b - a1*a4*log(TR) + a2*a3*t^b - a2*a4)*41.59 template. So I set

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

Because it does not have product between parameters, it seems to me much more natural to fit a model like this than the original model. It is more natural for me and easier on the computer, because it reduces the chance of numeric instabilities occurring.

In the end it looks 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

Please note that I used the same initial kicks as 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, has a value other than a1 in your parameterization. If there is a physical meaning to this parameter, it is good to see how to make the correct mathematical transformation to get the value of it in the context of your problem. The same goes for a2 , a3 and a4 .

    
12.05.2017 / 15:32