Fitting a physical model to a specific data using nls: over-parameterization or unidentifiable parameters?

Question

I have somewhat a complex physical model with five unknown parameters to fit, but no success so far.

I used nls2 first to get some estimates for the start values, but then nls, nlxb, and nlsLM all threw the famous "singular gradient error at initial parameter estimates" error.

For the start values for nls2, I extracted them from the literature, so I think that I have good starting values at least for nls2. The parameter estimates extracted from nls2 make quite sense physically as well; however, don't resolve the issue with the singular gradient matrix error.

Since it's a physical model, every coefficient has a physical meaning, and I prefer not to fix any of them.

I should also mention that all five unknown parameters in the model equation are positive and the shape parameter m can go up to 2.

Reading through many posts and trying different solution suggestions, I have come to conclusion that I have either over-parameterization or unidentifiable parameters problem.

My question is that should I stop trying to use nls with this specific model (with this many unknown parameters) or is there any way out?

I am quite new to topic, so any help, mathematically or code-wise, is greatly appreciated.

Here is my MWE:

# Data
x <-  c(0, 1000, 2000, 2500, 2750, 3000, 3250, 3500, 3750, 4000, 5000)
y <- c(1.0, 0.99, 0.98, 0.95, 0.795, 0.59, 0.35, 0.295, 0.175, 0.14, 0.095)

# Start values for nls2
bounds <- data.frame(a = c(0.8, 1.5), b = c(1e+5, 1e+7), c = c(0.4, 1.4), n = c(0.1, 2), m = c(0.1, 2))

# Model equation function
mod <- function(x, a, b, c, n, m){
  t <- b*85^n*exp(-c/0.0309)
  (1 - exp(-(a/(t*x))^m))
}

# # Model equation
# mod <- y ~ (1 - exp(-(a/(b*85^n*exp(-c/0.0309)*x))^m))

# Model fit with nls2 
fit2 <- nls2(y ~ mod(x, a, b, c, n, m), data = data.frame(x, y), start = bounds, algorithm = "brute-force")

# Model fit with nls
fit <- nls(y ~ mod(x, a, b, c, n, m), data = data.frame(x, y), start = coef(fit2))

user_190277 · Answer

The more I look at this the more confused I get, but I'm going to try again.

Looking again at your expression, we have the expression inside the exponential

-(a/(b*85^n*exp(-c/0.0309)*x))^m

We can rewrite this as


-( [a/(b*85^n*exp(-c/0.0309))] * 1/x )^m

(please check my algebra!) If this is correct, then that whole bold-faced blob doesn't affect the functional form of x — it all collapses to a single constant in the equation. (In other words, {a,b,c,n} are all jointly unidentifiable.) Lumping that stuff into a single parameter phi :

1 - exp(-(phi/x)^m)

phi is a shape parameter (has the same units as x, should be roughly the same magnitude as a typical value of x): let's try a starting value of 2500 (the mean value of x)
m is a shape parameter; we can't go too badly wrong starting from m==1

Now nls works fine without any extra help:

n1 <- nls(y~1 - exp(-(phi/x)^m), start=list(phi=2500,m=1), data=data.frame(x,y))

and gets phi=2935, m=6.49.

Plot predictions:

plot(x,y, ylim=c(0,1))
xvec <- seq(0, 5000, length=101)
lines(xvec, predict(n1, newdata=data.frame(x=xvec)))

Another way to think about what this curve is doing: we can transform the equation to -log(1-y) = phi^m*(1/x)^m: that is, -log(1-y) should follow a power-law curve with respect to 1/x.

Here's what that looks like:

plot(1/x, -log(1-y))
## curve() uses "x" as the current x-axis variable, i.e.
##  read "x" as "1/x" below.
with(as.list(coef(n1)), curve(phi^m*x^m, add=TRUE))

In this format, it appears to fit the central data well but fails for large values of 1/x (the x=0 point is missing here because it goes to infinity).

Fitting a physical model to a specific data using nls: over-parameterization or unidentifiable parameters?

1 Answers