[R] trying to obtain same nls parameters as in example

[Ben Bolker]

Pedro Mardones <mardones.p <at> gmail.com> writes:

> 
> Dear R-users;
> 
> I'm working with a a dataset that was previously used to fit a
> nonlinear model of the form:
> 
> Y ~ a * (1 + b * log(1 - c * X^d))
> 
> The parameters published elsewhere are:
> 
> a = 1.758863, b = .217217, c = .99031, and d = .054589
> 
> However, there is no way I can replicate this result. I've tried
> several options (including SAS) w/o success.
> 

  This is going to be a very difficult model to fit. 

## plot data and putative (old) fit
 oldparms <- list(a = 1.758863, b = .217217, c = .99031,  d = .054589)
plot(Y~X,data=test)
with(oldparms,curve(a * (1 + b * log(1 - c * x^d)),add=TRUE, col=2))

## compute RSS of old fit:

oldfit <- with(c(oldparms,as.list(test)),
               a * (1 + b * log(1 - c * X^d)))
oldresid <- test$Y-oldfit
sum(oldresid^2)  ## 0.263

 I used your nlsLM() code with a slight modification:

library(minpack.lm)
fit2 <- nlsLM(Y ~ a * (1 + b * log(1 - c * X^d)), data = test, 
  start = list(a = 1.75, b = .22, c = .99, d = .05),
    control = nls.lm.control(maxfev = 1000,maxiter=1000))

sum(residuals(fit2)^2)  ## 0.241
with(as.list(coef(fit2)),
     curve(a * (1 + b * log(1 - c * x^d)),add=TRUE, col=4))

## calculate correlation among parameters
cov2cor(vcov(fit2))

           a          b          c          d
a  1.0000000 -0.9999982  0.9999960 -0.9999991
b -0.9999982  1.0000000 -0.9999995  0.9999999
c  0.9999960 -0.9999995  1.0000000 -0.9999988
d -0.9999991  0.9999999 -0.9999988  1.0000000

  My conclusions: 
* nlsLM() finds a slightly better fit to the data than 
the quoted parameters;
* the parameters are all _extremely_ strongly
correlated with each other.  With a little bit more
algebra/calculus you could probably figure out why, by
Taylor expanding etc.: c is very close to 1, d is very close to zero ...