Hello colleagues,
I am attempting to determine the nonlinear least-squares estimates of the nonlinear model parameters using nls. I have come across a common problem that R users have reported when I attempt to fit a particular 3-parameter nonlinear function to my dataset:
Error in nls(r ~ tlm(a, N.fix, k, theta), data = tlm.data, start =
list(a = a.st, :
step factor 0.000488281 reduced below `minFactor' of 0.000976563
Despite modifying minFactor using nls.control, I am unable to counter the apparent singularity in the model fit. I have also tried changing the tolerance and start parameter values to no avail. If anyone can provide a relatively simple solution (perhaps adjusting the gradient, but I'm not sure how to do this), I would be most appreciative. My dataset is:
tlm.data
r N.fix
1 -0.52407085 76
2 0.10536052 45
3 -0.17435339 50
4 0.19415601 42
5 0.48701498 51
6 -0.50681760 83
7 -0.17435339 50
8 0.55278982 42
9 0.15219182 73
10 0.49899117 85
11 0.10821358 140
12 -0.83034830 156
13 -0.30748470 68
14 -0.22314355 50
15 0.04879016 40
16 -0.04879016 42
17 0.75377180 40
18 -0.12516314 85
19 -0.36624439 75
My function is:
tlm <- function(a,N,k,theta) (a*(1-((N/k)^theta)))
The nls fit I've coded is:
tlm.fit <- try(nls(r~tlm(a,N.fix,k,theta), data=tlm.data, start=list(a=a.st,k=k.st,theta=1),
trace=TRUE, control=nls.control(maxiter=6000,tol=1e-05,minFactor=1/1024)))
I'm using start values parsed in from another (previous, but not shown)
model fit. In this case,
a.st
[1] 0.3812922
k.st
[1] 64.66529
I happen to know the true values for the optimised parameters (from another application), but I can't get nls to reproduce them. They are:
a = 2.0466
k = 60.8275
theta = 0.2277
Any ideas?
Regards,
Corey Bradshaw
It is possible to fit this model to these data using nls as shown in the enclosed transcript. You have one conditionally linear parameter ('a') in the model so I used the plinear algorithm and I also generated analytic derivatives for the tls function using the deriv function.
There are several things to note:
- Your data are very noisy. It is not surprising that it is difficult to fit a 3-parameter nonlinear model to such data.
- The fitted model has negative values for both 'a' and 'theta'.
- The estimates are highly imprecise.
> summary(fm1)
Formula: r ~ tlm(N.fix, k, theta)
Parameters:
Estimate Std. Error t value Pr(>|t|)
k 49.0724 6.9153 7.096 2.53e-06
theta -4.6676 5.1805 -0.901 0.381
.lin -0.2333 0.1685 -1.385 0.185Residual standard error: 0.3662 on 16 degrees of freedom
Correlation of Parameter Estimates:
k theta
theta 0.8353
.lin -0.3458 -0.7104
R : Copyright 2004, The R Foundation for Statistical Computing Version 2.0.1 (2004-11-15), ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.
> dat <- read.table("/tmp/tlm.dat", header = TRUE)
> dat
r N.fix
1 -0.52407085 76
2 0.10536052 45
3 -0.17435339 50
4 0.19415601 42
5 0.48701498 51
6 -0.50681760 83
7 -0.17435339 50
8 0.55278982 42
9 0.15219182 73
10 0.49899117 85
11 0.10821358 140
12 -0.83034830 156
13 -0.30748470 68
14 -0.22314355 50
15 0.04879016 40
16 -0.04879016 42
17 0.75377180 40
18 -0.12516314 85
19 -0.36624439 75
> plot(r ~ N.fix, dat)
> tlm <- deriv( ~ 1-((N/k)^theta), c("k", "theta"), function(N, k, theta){})
> tlm
function (N, k, theta)
{
.expr1 <- N/k
.expr2 <- .expr1^theta
.value <- 1 - .expr2
.grad <- array(0, c(length(.value), 2), list(NULL, c("k",
"theta")))
.grad[, "k"] <- .expr1^(theta - 1) * (theta * (N/k^2))
.grad[, "theta"] <- -(.expr2 * log(.expr1))
attr(.value, "gradient") <- .grad
.value
}
> fm1 <- nls(r ~ tlm(N.fix, k, theta), dat, start=c(k = 65, theta = 1), alg =
> "plinear", trace = TRUE)
2.347334 : 65.0000000 1.0000000 0.3836488
2.233973 : 56.633165 -0.431863 -1.282247
2.187513 : 54.8608041 -1.3747423 -0.4908351
2.167521 : 52.8247351 -2.1635728 -0.3592099
2.155828 : 51.4475043 -2.8951236 -0.3003772
2.149315 : 50.4174840 -3.5536098 -0.2682023
2.146427 : 49.7548459 -4.0636580 -0.2500007
2.145514 : 49.3861622 -4.3828441 -0.2405601
2.145301 : 49.2043941 -4.5476404 -0.2362051
2.145260 : 49.1221764 -4.6225967 -0.2343549
2.145252 : 49.0870152 -4.6544862 -0.2335967
2.145251 : 49.0724421 -4.6676430 -0.2332897
> newdat <- list(N.fix = seq(40, 160))
> lines(newdat$N.fix, predict(fm1, newdat))
> q('no')
Process R finished at Wed Feb 23 11:20:23 2005
______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
