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
______________________________________________
R-help@r-project.org 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.
