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
???
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 [[alternative HTML version deleted]]
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
[[alternative HTML version deleted]]