padpadpadpad / nls.multstart Goto Github PK

Write vignette comparing fitting potential of nls and nlsLM with nls_multstart for multiple models.

does what it says on the tin.

Check model with more example fits.

Hey everyone, (first post here 😃 )

I have been running the same script once or twice already but this time, for some reason it does not work.
I have a very simple dataset with the speed of individuals at different temperature and I want to draw and fit a thermal response curve following the sharpe-Schoolfield model

`start_vals <- get_start_vals(speed$temp, speed$medianSpeed, model_name = 'sharpeschoolhigh_1981')

low_lims <- get_lower_lims(speed$temp, speed$medianSpeed, model_name = 'sharpeschoolhigh_1981')

upper_lims <- get_upper_lims(speed$temp, speed$medianSpeed, model_name = 'sharpeschoolhigh_1981')`

So this step works fine, but then I do:
fit <- nls_multstart(speed~sharpeschoolhigh_1981(temp = temperature, r_tref,e,eh,th, tref = 15), data = speed, iter = 500, start_lower = start_vals - 10, start_upper = start_vals + 10, lower = low_lims, upper = upper_lims, supp_errors = 'Y')

As simple as that. However, I still have the issue

Error in nls_multstart(speed ~ sharpeschoolhigh_1981(temp = temperature, : There must be as many parameter starting bounds as there are parameters

I tried to look at my start_vals vector and I got NA for 'eh'. I don't know if there is any issue related to that.

Any help would be more than welcome,

Thanks
Alann

I updated my tidyr the other day, and I get a warning every time I run nls_multstart.

Warning message:
All elements of `...` must be named.

This relates to the new syntax of the tidyr nest() command. I'll submit a PR to fix this.

Update: PR #16

I spotted your package when it popped up on my Twitter feed, and it looks really great, and it solves a problem which I've been looking to get to for my pharmacokinetic modelling package (https://github.com/mathesong/kinfitr). As such, I would very much like to make use of your package in mine. There are just a few changes/enhancements which I would like to make: I'd be happy to submit them as a pull request. Just wanted to check if they were ok for you first.

(by the way, I've never made a PR before, so please tell me if I should be doing anything differently here)

Data input argument

I would like to be able to include the data as vectors rather than as columns of a vector in some cases. I see that the minpack.lm package has achieved this by using parent.frame() as the default input argument. This allows that the data can be flexibly entered either in vector form in the formula, or in the data frame.

Parameter bounds

You implemented the parameter bounds input as a vector of lower1, upper1, lower2, upper2. This differs from the underlying minpack.lm input, and makes it a little bit harder to directly use the parameter bounds inputs which might have gone directly into minpack.lm as these inputs. Is there a special reason for putting them together? If not, would it be ok for me to separate the upper and lower input arguments for consistency with minpack.lm below?

runif shotgun

You implemented a shotgun approach selecting from a uniform distribution for choosing values for input arguments. Another possibility for scanning across the parameter space would be to divide each parameter's upper and lower starting parameter limit into x equally spaced units, and then fitting the model using each of these starting points. This would ensure that the model was fitted using starting parameters which covered the whole bounding space. I would like to include this option for having fewer iterations, but still covering the full space.

Again, the package looks great, and I'm really excited to try to implement it in my functions.

All the best,
Granville

I have been wandering around to find a good way for global fitting in R. I used nls.multstart and it works predictably (awesome package)
I want to use this package like I can use with nls in base R.

#Generating sample data
d <- transform(
    data.frame(
        x = seq(0, 1, len = 17),
        group = rep(c("A", "B", "B", "C"), len = 17)
    ),
    y = round(1 / (1.4 + x^ifelse(group == "A", 2.3, ifelse(group == "B", 3.1, 3.5))), 2),
    group = as.factor(group)
)

#Fitting data 
fit <- nls(y ~ 1 / (b + x^p[group]),
    data = d,
    start = list(b = 1, p = rep(3, length(levels(d$group))))
)

"The code is from here"
In the sample code above I can use "[]" to asign a parameter to each group while keeping "b"
as shared parameter.

I tried to follow the same pattern in nls.multstart but I get the following error

fitMulti <- nls.multstart::nls_multstart(
    data = d,
    y ~ 1 / (b + x^p[group]),
    iter = 200,
    start_lower = list(b = 0, p = rep(0, length(levels(d$group)))),
    start_upper = list(b = 5, p = rep(5, length(levels(d$group)))),
    supp_errors = "Y",
    convergence_count = 100,
    na.action = na.omit
)

#ouput
 Error in `as_tibble()`:
! Column names `b` and `p` must not be duplicated.
Use `.name_repair` to specify repair.
Caused by error in `repaired_names()`:
! Names must be unique.
✖ These names are duplicated:
  * "b" at locations 1 and 3.
  * "p" at locations 2 and 4.
Run `rlang::last_trace()` to see where the error occurred.

Is there a way to use nls.multstart to do global fitting in this context?

Hi Daniel & Granville
@mathesong
@padpadpadpad

I have these issues:

• Every time I call nls_multstart, I get different results for the fitting parameters on the same exact data. How can I get this fixed and get all time similar results for the fitting parameters ( eta , ntotperl , sig ).

• As you see from the image below, I am fitting two data (blueData, redData) in for loop and get the fitting parameters for each separately. Expecting I get similar fitting parameters for (eta , ntotperl , sig) for both data. As you see in the image, those two curves very similar except they are pi phase shifted from one another.

• I noticed the fitting sometimes takes the value I assigned in 'lower' in nls_multstart function as shown in the image below , ntotperl=0.00e+00
  If I don't specify values in "lower", I get negative valuse for ( eta , ntotperl , sig) which are not realistic.
  For both data (blueData, redData), the expected positive values for (eta , ntotperl , sig) are eta =~ 2.40e5, ntotperl =~ 2.34e-9, sig = ~8.0e10 and they could be different but they need to be consistent every time I run nls_multstart
  For beta1 and beta2 they are random phase or angels in radiance and they can take any values. p1 is the max-min for each data ~ 0.002 -0.004

Here are the data, function , and nls_multstart input parametrs
your help is really appreciated

Thank you
Abdo

sitedata

OPT_blueData:
0.0029621019 0.0031441264 0.0033358000 0.0034986460 0.0036702868 0.0038132915 0.0039111097 0.0040373835 0.0040803743 0.0040709897 0.0040105137 0.0038815037 0.0036576319 0.0034387452 0.0031434025 0.0028047877 0.0023944182 0.0019696996 0.0015332042 0.0011084089 0.0007231027 0.0003991168 0.0001695899 0.0000978363
OPT_redData:
1.298375e-03 1.122535e-03 9.253372e-04 7.426770e-04 5.618236e-04 3.994846e-04 2.674853e-04 1.632676e-04 1.043999e-04 9.799412e-05 1.429881e-04 2.435567e-04 4.067242e-04 6.214411e-04 8.943345e-04 1.213109e-03 1.571810e-03 1.957492e-03 2.328627e-03 2.716439e-03 3.082478e-03 3.345800e-03 3.505097e-03 3.535086e-03
Volt:
-4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0
lam:
1264

custom function

`
mzmfunc <- function(p1, beta1, beta2, eta, ntotperl , sig , Volt1, Volt2,lamb) {

lamb<- lamb*1e-9

l<-0.63e-2

co_mod<- 1

capperl<-0.26e-9

cap <- capperllco_mod

vbi_va <- 0.4

m_va <- 0.3

mgrad <- (1.0/m_va) - 2.0

qo <- -1*((mgrad+2)/(mgrad+1))*exp(log(vbi_va)/(mgrad+2))*cap

qt <- qo*exp(((mgrad+1)/(mgrad+2))*log(vbi_va))

ntot <- ntotperllco_mod*co_mod

qcap1<- qt*exp(((mgrad+1)/(mgrad+2))*log(1- Volt1/vbi_va))

phi1<- (-2*pi/lamb)etaqcap1

qcap2<- qt*exp(((mgrad+1)/(mgrad+2))*log(1- Volt2/vbi_va))

phi2<- (-2*pi/lamb)etaqcap2

A1<- sqrt(abs(p1*( 1-sig*(ntot+qcap1))))

A2<- sqrt(abs(p1*( 1-sig*(ntot+qcap2))))

Eo1<- -1i * 0.5 ( A1exp(-1i * (phi1+beta1)) + A2*exp(-1i * (phi2+beta2)) )

Eo2<- 0.5*( A1*exp(-1i (phi1+beta1)) - A2exp(-1i * (phi2+beta2)) )

pow1<- abs(Eo1*Conj(Eo1))

pow2<- abs(Eo2*Conj(Eo2))

return(pow1)
}
`

nls_multstart

`
result1 <- nls_multstart(OPT~ mzmfunc( p1, beta1, beta2, eta, ntotperl, sig , Volt1 = 0, Volt2 = Volt, lamb =lam),

                       data = sitedata,

                       iter = 1000,

                       start_lower = c(p1 = 0.001, beta1 = 0, beta2 = 0, eta = 0 , ntotperl = 0, sig = 0),

                       start_upper = c(p1 = 0.006 , beta1 = 2*pi, beta2 = 2*pi, eta = 3e5 , ntotperl = 1e-9, sig =10e10),

                       supp_errors = 'Y',

                       convergence_count = 300,

                       na.action = na.omit,

                       lower = c(p1 = 0.0001,beta1=0, beta2 = 0, eta =0, ntotperl = 0, sig = 0))

`

After the implementation of the updates suggested by @mathesong, I think the package should get submitted to CRAN. Just a single function it does a valuable job for the R community.

If only so that the package can start earning citation when me and my lab use it in the methods of our manuscripts! Free references are not to be missed. And it will also help publicise the package so people actually know about it!

Links for preparing a package for CRAN are from Karl Broman, Hadley Wickham and CRAN

Hello,
First I would like to thank you for this great work .
I am trying to use nls.multstart to fit sine like curve, but I getting error that "There must be as many parameter starting bounds as there are parameters". I double checked all parameters and they are same as the bounds parameters. Below is my function and code I used:

As you can see below, I am trying to fit OPT vs. Volts. The fitting parameters are six: p1, beta1, beta2, eta, ntot , sig
The actual data has multiple sets/sites of similar data shown below, that is why I want to use your function which supports multiple starting values.
Is it too many parameters to fit and that is why it fails? Is it possible to make this fitting happen ? Could you please help me to make this working ?

Thank you
Abdo

function to fit

mzmfunc <- function(p1, beta1, beta2, eta, ntot , sig , Volt1, Volt2,lamb) {

l<-0.63e-2
co_mod<- 1
capperl<-0.26e-9
cap <- capperl*l*co_mod
vbi_va <- 0.4
m_va <- 0.3
mgrad <- (1.0/m_va) - 2.0
qo <- -1*((mgrad+2)/(mgrad+1))*exp(log(vbi_va)/(mgrad+2))*cap
qt <- qo*exp(((mgrad+1)/(mgrad+2))*log(vbi_va))
ntotperl <- 2.34e-9
ntoti <- ntotperl*l*co_mod*co_mod
qcap1<- qt*exp(((mgrad+1)/(mgrad+2))*log(1- Volt1/vbi_va))
phi1<- (-2*pi/lamb)*eta*qcap1
qcap2<- qt*exp(((mgrad+1)/(mgrad+2))*log(1- Volt2/vbi_va))
phi2<- (-2*pi/lamb)*eta*qcap2

A1<- sqrt(abs(p1*( 1-sig*(ntot+qcap1))))
A2<- sqrt(abs(p1*( 1-sig*(ntot+qcap2))))
Eo1<- -1i*0.5*( A1*exp(-1i*(phi1+beta1)) + A2*exp(-1i*(phi2+beta2)) )
Eo2<- 0.5*( A1*exp(-1i*(phi1+beta1)) - A2*exp(-1i*(phi2+beta2)) )
pow1<- abs(Eo1*Conj(Eo1))
pow2<- abs(Eo2*Conj(Eo2))
return(pow1)

}

using nls_multstart

result1 <- nls_multstart(OPT ~ mzmfunc(p1, beta1, eta, ntot , sig , Volt1=0, Volt2=Volt, lamb=1.262e-6),
data = sitedata,
iter = 250,
start_lower = c(p1= 0.01, beta1 = 0 , beta2 = 0, eta=1.4e5 , ntot= 1.4742e-11, sig= 1e10 ),
start_upper = c(p1= 0.04, beta1 = 4 , beta2 = 4, eta=10.4e5 , ntot= 10.4742e-11, sig= 10e10 ),
supp_errors = 'Y',
convergence_count = 200,
na.action = na.omit,
lower = c( p1=0.01, beta1 = 0, beta2 = 0, eta=2.4e5, ntot= 1.4742e-11, sig= 1e10))

sitedata
Volt
-4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0
OPT
0.005701643 0.006880183 0.008212966 0.009582967 0.011132687 0.012385114 0.013740420 0.014845677 0.016180800 0.017108030 0.017901936 0.018365383 0.018365383 0.018395009 0.017926686 0.016900517 0.015848932 0.014197112 0.012705741 0.010514775 0.008077928

Hey Dan,

I tried out the package on my specific problem, and it seems to work really well! However, passing arguments through to minpack.lm turned out to be a little bit trickier than expected. Upper and lower worked a charm, but weights gave me errors each time. As far as I can tell, this is because minpack.lm underneath is doing some kind of non-standard evaluation on weights, and it doesn't work when it's obscured a level by nls.multstart.

I just tried out a very simple solution: adding weights as an optional argument to nls.multstart fixes it completely. If not specified, it works as usual. If specified, then it allows one to do weighted nls.

So I'll quickly submit a PR with that.