Error in chol.default

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Error in chol.default

Luigi
Hello,
I am trying to fit a Richards model to some cumulative incidence data
of infection. I got this example:
```
rich = function(p, x) {
  a = p["curvature"]
  k = p["finalPop"]
  r = p["growthRate"]
  y = r * x * (1-(x/k)^a)
  return(y)
}
ricky = function(p, x, y) p$r * x * (1-(x/p$k)^p$a) -y
# values
Y <- c(41,   41,   41,   41,   41,   41,   45,   62,  121,  198,  275,
 288,  859, 1118)
X <- 1:14
a = 1
k = 83347
r = 40
fit = rich(c(curvature=a, finalPop=k, growthRate=r), X)
plot(Y ~ X,
     col = "red", lwd = 2,
     main = "Richards model",
     xlab = expression(bold("Days")),
     ylab = expression(bold("Cases")))
points(X, fit, type = "l", lty = 2, lwd = 2)
library("minpack.lm")
o <- nls.lm(par = list(a=a, k=k, r=r), fn = ricky, x = X, y = Y)
summary(o)
```
When I run the optimization (given that I can't find parameters that
fit the data by eyeball), I get the error:
```
Error in chol.default(object$hessian) :
  the leading minor of order 1 is not positive definite
```
What is this about?
Thank you




--
Best regards,
Luigi

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: Error in chol.default

J C Nash
The error msg says it all if you know how to read it.

> When I run the optimization (given that I can't find parameters that
> fit the data by eyeball), I get the error:
> ```
> Error in chol.default(object$hessian) :
>   the leading minor of order 1 is not positive definite


Your Jacobian (derivatives of the residual function w.r.t. the parameters)
is singular -- rather spectacularly so.

Try the short addition to your code to use analytic derivatives:

rich = function(p, x) {
  a = p["curvature"]
  k = p["finalPop"]
  r = p["growthRate"]
  y = r * x * (1-(x/k)^a)
  return(y)
}
ricky = function(p, x, y) p$r * x * (1-(x/p$k)^p$a) -y
# values
Y <- c(41,   41,   41,   41,   41,   41,   45,   62,  121,  198,  275,
       288,  859, 1118)
X <- 1:14
a = 1
k = 83347
r = 40
fit = rich(c(curvature=a, finalPop=k, growthRate=r), X)
plot(Y ~ X,
     col = "red", lwd = 2,
     main = "Richards model",
     xlab = expression(bold("Days")),
     ylab = expression(bold("Cases")))
points(X, fit, type = "l", lty = 2, lwd = 2)
library("minpack.lm")
o <- nls.lm(par = list(a=a, k=k, r=r), fn = ricky, x = X, y = Y)
summary(o)
print(o1)
library(nlsr)
xy=data.frame(x=X, y=Y)
rcky2 = y ~ r * x * (1-(x/k)^a) -y
o1 <- nlxb(rcky2, start = c(a=a, k=k, r=r), data=xy, trace=TRUE)


You should find

> o1
nlsr object: x
residual sumsquares =  3161769  on  14 observations
    after  8    Jacobian and  13 function evaluations
  name            coeff          SE       tstat      pval      gradient    JSingval
a                294.113            NA         NA         NA           0       31.86
k                  83347            NA         NA         NA           0           0
r                74.9576            NA         NA         NA  -6.064e-05           0
>
Note the singular values of the Jacobian. Actual zeros!

Even so, nlsr had managed to make some progress.

nls() and nls.lm() use approximate derivatives. Often that's fine (and it is more flexible
in handling awkward functions), but a lot of the time it is NOT OK.

JN

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.