# [R] nls with plinear and function on RHS

Keith Jewell k.jewell at campden.co.uk
Thu Oct 2 14:03:39 CEST 2008

```Dear R gurus,

As part of finding initial values for a much more complicated fit I want to
fit a function of the form y ~ a + bx + cx^d to fairly "noisy" data and have
hit some problems.

To demonstrate the specific R-related problem, here is an idealised data
set, smaller and better fitting than reality:
# idealised data set
aDF <- data.frame( x= c(1.80, 9.27, 6.48, 2.61, 9.86, 5.93, 6.76, 5.52,
6.06, 8.62),
y= c(24.77, 2775.07, 895.15, 60.73, 3373.57, 677.82, 1021.92, 542.84,
725.25, 2200.04))

And here are some starting values, far better than I'd have in reality
# good starting values
startL <- list(Ca=4, Cb=3, Cc=3, Cd=3)

In these idealised circumstances nls converges using the default algorithm.
# nls, default algorithm
nls(y ~ Ca + Cb*x + Cc*x^Cd, data=aDF, start=startL)

Unfortunately, in reality it often fails to converge. This model is linear
in a, b and c so I've used the plinear algorithm.
# nls, plinear algorithm, explicit RHS
nls(y ~ cbind(Ca=1,Cb=x, Cc=x^Cd), data=aDF, start=startL["Cd"],
algorithm="plinear")

This converges much more often but sometimes crashes with the error message
"Error in numericDeriv(form[[3]], names(ind), env) :
Missing value or an infinity produced when evaluating the model"

I deduce that it is failing in the numerical differentiation of x^Cd (but
don't know why), so I thought I'd avoid the numerical differentiation by
putting the RHS in a function to which I could (later) add a 'gradient'
attribute
# function to provide RHS
bFunc <- function(x, Ca, Cb, Cc, Cd) cbind(Ca=1,Cb=x, Cc=x^Cd)

# nls, plinear algorithm, RHS from function
nls(y ~ bFunc(x, Ca, Cb, Cc, Cd), data=aDF, start=startL["Cd"],
algorithm="plinear")

However, this gives me
"Error in nls(y ~ bFunc(x, Ca, Cb, Cc, Cd), data = aDF, start =
startL["Cd"],  :
parameters without starting value in 'data': Ca, Cb, Cc"

Can anyone tell me
a) why putting the RHS into a function "broke" the plinear algorithm
b) if there's a better approach to my problem

Keith Jewell
-----------------
I'm using V2.7.2...
> sessionInfo()
R version 2.7.2 (2008-08-25)
i386-pc-mingw32

locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United
Kingdom.1252;LC_MONETARY=English_United
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

attached base packages:
[1] stats     graphics  grDevices datasets  tcltk     utils     methods
base

other attached packages: