[R] How to coerce a parameter in nls?
Jianling Fan
fanjianling at gmail.com
Tue Sep 22 17:12:13 CEST 2015
Hello Gabor,
It is very kind of you to reply and give suggestion so rapid. I will
try to learn and use it.
Thanks very much for your help!
Best regards,
Jianling
On 22 September 2015 at 06:45, Gabor Grothendieck
<ggrothendieck at gmail.com> wrote:
> Or if you really can't bear to write out 20 terms have R do it for you:
>
> # number of terms is the number of unique values in ref column
> nterms <- length(unique(dproot$ref))
>
> dproot2 <- do.call(data.frame, transform(dproot, ref = outer(dproot$ref,
> seq(nterms), "==") + 0))
>
> # construct the formula as a string
> terms <- paste( sprintf("Rm%d*ref.%d", 1:nterms, 1:nterms), collapse = "+")
> fo <- sprintf("den ~ (%s)/(1+(depth/d50)^c)", terms)
>
> library(nlmrt)
> fm <- nlxb(fo, data = dproot2, masked = "Rm6",
> start = c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65, Rm5=1.01, Rm6=1,
> d50=20, c=-1))
>
>
> On Tue, Sep 22, 2015 at 7:04 AM, Gabor Grothendieck
> <ggrothendieck at gmail.com> wrote:
>>
>> Just write out the 20 terms.
>>
>> On Mon, Sep 21, 2015 at 10:26 PM, Jianling Fan <fanjianling at gmail.com>
>> wrote:
>>>
>>> Hello, Gabor,
>>>
>>> Thanks again for your suggestion. And now I am trying to improve the
>>> code by adding a function to replace the express "Rm1 * ref.1 + Rm2 *
>>> ref.2 + Rm3 * ref.3 + Rm4 * ref.4 + Rm5 * ref.5 + Rm6 * ref.6" because
>>> I have some other dataset need to fitted to the same model but with
>>> more groups (>20).
>>>
>>> I tried to add the function as:
>>>
>>> denfun<-function(i){
>>> for(i in 1:6){
>>> Rm<-sum(Rm[i]*ref.i)
>>> return(Rm)}
>>> }
>>>
>>> but I got another error when I incorporate this function into my
>>> regression:
>>>
>>> >fitdp1<-nlxb(den ~ denfun(6)/(1+(depth/d50)^c),
>>> data = dproot2,
>>> start = c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65,
>>> Rm5=1.01, Rm6=1, d50=20, c=-1),
>>> masked = "Rm6")
>>>
>>> Error in deriv.default(parse(text = resexp), names(start)) :
>>> Function 'denfun' is not in the derivatives table
>>>
>>> I think there must be something wrong with my function. I tried some
>>> times but am not sure how to improve it because I am quite new to R.
>>>
>>> Could anyone please give me some suggestion.
>>>
>>> Thanks a lot!
>>>
>>>
>>> Jianling
>>>
>>>
>>> On 22 September 2015 at 00:43, Gabor Grothendieck
>>> <ggrothendieck at gmail.com> wrote:
>>> > Express the formula in terms of simple operations like this:
>>> >
>>> > # add 0/1 columns ref.1, ref.2, ..., ref.6
>>> > dproot2 <- do.call(data.frame, transform(dproot, ref =
>>> > outer(dproot$ref,
>>> > seq(6), "==") + 0))
>>> >
>>> > # now express the formula in terms of the new columns
>>> > library(nlmrt)
>>> > fitdp1<-nlxb(den ~ (Rm1 * ref.1 + Rm2 * ref.2 + Rm3 * ref.3 + Rm4 *
>>> > ref.4 +
>>> > Rm5 * ref.5 + Rm6 * ref.6)/(1+(depth/d50)^c),
>>> > data = dproot2,
>>> > start = c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65, Rm5=1.01,
>>> > Rm6=1,
>>> > d50=20, c=-1),
>>> > masked = "Rm6")
>>> >
>>> > where we used this input:
>>> >
>>> > Lines <- " depth den ref
>>> > 1 20 0.5730000 1
>>> > 2 40 0.7800000 1
>>> > 3 60 0.9470000 1
>>> > 4 80 0.9900000 1
>>> > 5 100 1.0000000 1
>>> > 6 10 0.6000000 2
>>> > 7 20 0.8200000 2
>>> > 8 30 0.9300000 2
>>> > 9 40 1.0000000 2
>>> > 10 20 0.4800000 3
>>> > 11 40 0.7340000 3
>>> > 12 60 0.9610000 3
>>> > 13 80 0.9980000 3
>>> > 14 100 1.0000000 3
>>> > 15 20 3.2083491 4
>>> > 16 40 4.9683383 4
>>> > 17 60 6.2381133 4
>>> > 18 80 6.5322348 4
>>> > 19 100 6.5780660 4
>>> > 20 120 6.6032064 4
>>> > 21 20 0.6140000 5
>>> > 22 40 0.8270000 5
>>> > 23 60 0.9500000 5
>>> > 24 80 0.9950000 5
>>> > 25 100 1.0000000 5
>>> > 26 20 0.4345774 6
>>> > 27 40 0.6654726 6
>>> > 28 60 0.8480684 6
>>> > 29 80 0.9268951 6
>>> > 30 100 0.9723207 6
>>> > 31 120 0.9939966 6
>>> > 32 140 0.9992400 6"
>>> >
>>> > dproot <- read.table(text = Lines, header = TRUE)
>>> >
>>> >
>>> >
>>> > On Mon, Sep 21, 2015 at 12:22 PM, Jianling Fan <fanjianling at gmail.com>
>>> > wrote:
>>> >>
>>> >> Thanks Prof. Nash,
>>> >>
>>> >> Sorry for late reply. I am learning and trying to use your nlmrt
>>> >> package since I got your email. It works good to mask a parameter in
>>> >> regression but seems does work for my equation. I think the problem is
>>> >> that the parameter I want to mask is a group-specific parameter and I
>>> >> have a "[]" syntax in my equation. However, I don't have your 2014
>>> >> book on hand and couldn't find it in our library. So I am wondering if
>>> >> nlxb works for group data?
>>> >> Thanks a lot!
>>> >>
>>> >> following is my code and I got a error form it.
>>> >>
>>> >> > fitdp1<-nlxb(den~Rm[ref]/(1+(depth/d50)^c),data=dproot,
>>> >> + start =c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65,
>>> >> Rm5=1.01, Rm6=1, d50=20, c=-1),
>>> >> + masked=c("Rm6"))
>>> >>
>>> >> Error in deriv.default(parse(text = resexp), names(start)) :
>>> >> Function '`[`' is not in the derivatives table
>>> >>
>>> >>
>>> >> Best regards,
>>> >>
>>> >> Jianling
>>> >>
>>> >>
>>> >> On 20 September 2015 at 12:56, ProfJCNash <profjcnash at gmail.com>
>>> >> wrote:
>>> >> > I posted a suggestion to use nlmrt package (function nlxb to be
>>> >> > precise),
>>> >> > which has masked (fixed) parameters. Examples in my 2014 book on
>>> >> > Nonlinear
>>> >> > parameter optimization with R tools. However, I'm travelling just
>>> >> > now,
>>> >> > or
>>> >> > would consider giving this a try.
>>> >> >
>>> >> > JN
>>> >> >
>>> >> >
>>> >> > On 15-09-20 01:19 PM, Jianling Fan wrote:
>>> >> >>
>>> >> >> no, I am doing a regression with 6 group data with 2 shared
>>> >> >> parameters
>>> >> >> and 1 different parameter for each group data. the parameter I want
>>> >> >> to
>>> >> >> coerce is for one group. I don't know how to do it. Any suggestion?
>>> >> >>
>>> >> >> Thanks!
>>> >> >>
>>> >> >> On 19 September 2015 at 13:33, Jeff Newmiller
>>> >> >> <jdnewmil at dcn.davis.ca.us>
>>> >> >> wrote:
>>> >> >>>
>>> >> >>> Why not rewrite the function so that value is not a parameter?
>>> >> >>>
>>> >> >>>
>>> >> >>>
>>> >> >>> ---------------------------------------------------------------------------
>>> >> >>> Jeff Newmiller The ..... .....
>>> >> >>> Go
>>> >> >>> Live...
>>> >> >>> DCN:<jdnewmil at dcn.davis.ca.us> Basics: ##.#. ##.#.
>>> >> >>> Live
>>> >> >>> Go...
>>> >> >>> Live: OO#.. Dead: OO#..
>>> >> >>> Playing
>>> >> >>> Research Engineer (Solar/Batteries O.O#. #.O#.
>>> >> >>> with
>>> >> >>> /Software/Embedded Controllers) .OO#. .OO#.
>>> >> >>> rocks...1k
>>> >> >>>
>>> >> >>>
>>> >> >>>
>>> >> >>> ---------------------------------------------------------------------------
>>> >> >>> Sent from my phone. Please excuse my brevity.
>>> >> >>>
>>> >> >>> On September 18, 2015 9:54:54 PM PDT, Jianling Fan
>>> >> >>> <fanjianling at gmail.com> wrote:
>>> >> >>>>
>>> >> >>>> Hello, everyone,
>>> >> >>>>
>>> >> >>>> I am using a nls regression with 6 groups data. I am trying to
>>> >> >>>> coerce
>>> >> >>>> a parameter to 1 by using a upper and lower statement. but I
>>> >> >>>> always
>>> >> >>>> get an error like below:
>>> >> >>>>
>>> >> >>>> Error in ifelse(internalPars < upper, 1, -1) :
>>> >> >>>> (list) object cannot be coerced to type 'double'
>>> >> >>>>
>>> >> >>>> does anyone know how to fix it?
>>> >> >>>>
>>> >> >>>> thanks in advance!
>>> >> >>>>
>>> >> >>>> My code is below:
>>> >> >>>>
>>> >> >>>>
>>> >> >>>>
>>> >> >>>>> dproot
>>> >> >>>>
>>> >> >>>> depth den ref
>>> >> >>>> 1 20 0.5730000 1
>>> >> >>>> 2 40 0.7800000 1
>>> >> >>>> 3 60 0.9470000 1
>>> >> >>>> 4 80 0.9900000 1
>>> >> >>>> 5 100 1.0000000 1
>>> >> >>>> 6 10 0.6000000 2
>>> >> >>>> 7 20 0.8200000 2
>>> >> >>>> 8 30 0.9300000 2
>>> >> >>>> 9 40 1.0000000 2
>>> >> >>>> 10 20 0.4800000 3
>>> >> >>>> 11 40 0.7340000 3
>>> >> >>>> 12 60 0.9610000 3
>>> >> >>>> 13 80 0.9980000 3
>>> >> >>>> 14 100 1.0000000 3
>>> >> >>>> 15 20 3.2083491 4
>>> >> >>>> 16 40 4.9683383 4
>>> >> >>>> 17 60 6.2381133 4
>>> >> >>>> 18 80 6.5322348 4
>>> >> >>>> 19 100 6.5780660 4
>>> >> >>>> 20 120 6.6032064 4
>>> >> >>>> 21 20 0.6140000 5
>>> >> >>>> 22 40 0.8270000 5
>>> >> >>>> 23 60 0.9500000 5
>>> >> >>>> 24 80 0.9950000 5
>>> >> >>>> 25 100 1.0000000 5
>>> >> >>>> 26 20 0.4345774 6
>>> >> >>>> 27 40 0.6654726 6
>>> >> >>>> 28 60 0.8480684 6
>>> >> >>>> 29 80 0.9268951 6
>>> >> >>>> 30 100 0.9723207 6
>>> >> >>>> 31 120 0.9939966 6
>>> >> >>>> 32 140 0.9992400 6
>>> >> >>>>
>>> >> >>>>> fitdp<-nls(den~Rm[ref]/(1+(depth/d50)^c),data=dproot,
>>> >> >>>>
>>> >> >>>> + start = list(Rm=c(1.01, 1.01, 1.01, 6.65,1.01,1), d50=20,
>>> >> >>>> c=-1))
>>> >> >>>>>
>>> >> >>>>> summary(fitdp)
>>> >> >>>>
>>> >> >>>>
>>> >> >>>> Formula: den ~ Rm[ref]/(1 + (depth/d50)^c)
>>> >> >>>>
>>> >> >>>> Parameters:
>>> >> >>>> Estimate Std. Error t value Pr(>|t|)
>>> >> >>>> Rm1 1.12560 0.07156 15.73 3.84e-14 ***
>>> >> >>>> Rm2 1.57643 0.11722 13.45 1.14e-12 ***
>>> >> >>>> Rm3 1.10697 0.07130 15.53 5.11e-14 ***
>>> >> >>>> Rm4 7.23925 0.20788 34.83 < 2e-16 ***
>>> >> >>>> Rm5 1.14516 0.07184 15.94 2.87e-14 ***
>>> >> >>>> Rm6 1.03658 0.05664 18.30 1.33e-15 ***
>>> >> >>>> d50 22.69426 1.03855 21.85 < 2e-16 ***
>>> >> >>>> c -1.59796 0.15589 -10.25 3.02e-10 ***
>>> >> >>>> ---
>>> >> >>>> Signif. codes: 0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1
>>> >> >>>>
>>> >> >>>> Residual standard error: 0.1094 on 24 degrees of freedom
>>> >> >>>>
>>> >> >>>> Number of iterations to convergence: 8
>>> >> >>>> Achieved convergence tolerance: 9.374e-06
>>> >> >>>>
>>> >> >>>>> fitdp1<-nls(den~Rm[ref]/(1+(depth/d50)^c),data=dproot,
>>> >> >>>>
>>> >> >>>> algorithm="port",
>>> >> >>>> + start = list(Rm=c(1.01, 1.01, 1.01, 6.65, 1.01, 1), d50=20,
>>> >> >>>> c=-1),
>>> >> >>>> + lower = list(Rm=c(1.01, 1.01, 1.01, 6.65, 1.01, 1), d50=20,
>>> >> >>>> c=-1),
>>> >> >>>> + upper = list(Rm=c(2.1, 2.2, 2.12, 12.5, 2.3, 1), d50=50, c=1))
>>> >> >>>>
>>> >> >>>> Error in ifelse(internalPars < upper, 1, -1) :
>>> >> >>>> (list) object cannot be coerced to type 'double'
>>> >> >>>>
>>> >> >>>> ______________________________________________
>>> >> >>>> R-help at 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.
>>> >> >>>
>>> >> >>>
>>> >> >>
>>> >> >>
>>> >> >>
>>> >> >
>>> >> > ______________________________________________
>>> >> > R-help at 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.
>>> >>
>>> >>
>>> >>
>>> >> --
>>> >> Jianling Fan
>>> >> 樊建凌
>>> >>
>>> >> ______________________________________________
>>> >> R-help at 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.
>>> >
>>> >
>>> >
>>> >
>>> > --
>>> > Statistics & Software Consulting
>>> > GKX Group, GKX Associates Inc.
>>> > tel: 1-877-GKX-GROUP
>>> > email: ggrothendieck at gmail.com
>>>
>>>
>>>
>>> --
>>> Jianling Fan
>>> 樊建凌
>>
>>
>>
>>
>> --
>> Statistics & Software Consulting
>> GKX Group, GKX Associates Inc.
>> tel: 1-877-GKX-GROUP
>> email: ggrothendieck at gmail.com
>
>
>
>
> --
> Statistics & Software Consulting
> GKX Group, GKX Associates Inc.
> tel: 1-877-GKX-GROUP
> email: ggrothendieck at gmail.com
--
Jianling Fan
樊建凌
More information about the R-help
mailing list