nls() syntax

5 messages · Rolf Turner, Gabor Grothendieck, Duncan Murdoch +1 more

Original

1

5

Rolf Turner

Fri, Dec 11, 2020 3:50 PM #

I want to fit a model y = x/(x-a) where the value of a depends
on the level of a factor z.  I cannot figure out an appropriate
syntax for nls().  The "parameter" a (to be estimated) should be a
vector of length equal to the number of levels of z.

I tried:

strt <- rep(3,length(levels(z))
names(strt=levels(z)
fit <- nls(y ~ x/(x - a[z]),start=strt,data=xxx)

but of course got an error:

I keep thinking that there is something obvious that I should
be doing, but I can't work out what it is.

Does there *exist* an appropriate syntax for doing what I want
to do?  Can anyone enlighten me?  The data set "xxx" is given
in dput() form at the end of this message.

cheers,

Rolf Turner

Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

Data set "xxx":

structure(list(x = c(30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 
130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250, 
30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 
170, 180, 190, 200, 210, 220, 230, 240, 250, 30, 40, 50, 60, 
70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 
200, 210, 220, 230, 240, 250, 30, 40, 50, 60, 70, 80, 90, 100, 
110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 
240, 250, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 
150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250), y = c(1.27, 
1.16, 1.19, 1.15, 1.09, 1.07, 1.07, 1.05, 1.07, 1.03, 1.05, 1.07, 
1.06, 1.03, 1.05, 1.04, 1.03, 1.03, 1.03, 1.02, 1.02, 1.01, 1.01, 
1.21, 1.15, 1.1, 1.1, 1.06, 1.06, 1.05, 1.03, 1.07, 1.04, 1.04, 
1.02, 1.04, 1.02, 1.04, 1.03, 1.01, 1.03, 1.01, 1, 1.02, 1.03, 
1.02, 1.42, 1.27, 1.23, 1.14, 1.17, 1.08, 1.11, 1.06, 1.07, 1.08, 
1.06, 1.07, 1.04, 1.03, 1.07, 1.04, 1.03, 1.03, 1.03, 1.04, 1.03, 
1.03, 1.04, 1.85, 1.41, 1.35, 1.21, 1.22, 1.15, 1.14, 1.07, 1.1, 
1.09, 1.1, 1.09, 1.08, 1.08, 1.09, 1.09, 1.07, 1.06, 1.03, 1.08, 
1.05, 1.02, 1.05, 1.99, 1.6, 1.44, 1.4, 1.24, 1.3, 1.21, 1.23, 
1.18, 1.18, 1.12, 1.15, 1.09, 1.07, 1.13, 1.1, 1.05, 1.13, 1.09, 
1.03, 1.11, 1.07, 1.05), z = structure(c(1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), .Label = c("p1", 
"p2", "p3", "p4", "p5"), class = "factor")), class = "data.frame", row.names = c(NA, 
-115L))

Gabor Grothendieck

Fri, Dec 11, 2020 4:20 PM #

The start= argument should be as follows:

 nls(y ~ x/(x - a[z]),start=list(a = strt),data=xxx)

On Fri, Dec 11, 2020 at 6:51 PM Rolf Turner <r.turner at auckland.ac.nz> wrote:



I want to fit a model y = x/(x-a) where the value of a depends
on the level of a factor z.  I cannot figure out an appropriate
syntax for nls().  The "parameter" a (to be estimated) should be a
vector of length equal to the number of levels of z.

I tried:

strt <- rep(3,length(levels(z))
names(strt=levels(z)
fit <- nls(y ~ x/(x - a[z]),start=strt,data=xxx)

but of course got an error:

Error in nls(y ~ x/(x - a[z]), start = strt, data = xxx) :
  parameters without starting value in 'data': a

I keep thinking that there is something obvious that I should
be doing, but I can't work out what it is.

Does there *exist* an appropriate syntax for doing what I want
to do?  Can anyone enlighten me?  The data set "xxx" is given
in dput() form at the end of this message.

cheers,

Rolf Turner

--
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

Data set "xxx":

structure(list(x = c(30, 40, 50, 60, 70, 80, 90, 100, 110, 120,
130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250,
30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160,
170, 180, 190, 200, 210, 220, 230, 240, 250, 30, 40, 50, 60,
70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190,
200, 210, 220, 230, 240, 250, 30, 40, 50, 60, 70, 80, 90, 100,
110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230,
240, 250, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140,
150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250), y = c(1.27,
1.16, 1.19, 1.15, 1.09, 1.07, 1.07, 1.05, 1.07, 1.03, 1.05, 1.07,
1.06, 1.03, 1.05, 1.04, 1.03, 1.03, 1.03, 1.02, 1.02, 1.01, 1.01,
1.21, 1.15, 1.1, 1.1, 1.06, 1.06, 1.05, 1.03, 1.07, 1.04, 1.04,
1.02, 1.04, 1.02, 1.04, 1.03, 1.01, 1.03, 1.01, 1, 1.02, 1.03,
1.02, 1.42, 1.27, 1.23, 1.14, 1.17, 1.08, 1.11, 1.06, 1.07, 1.08,
1.06, 1.07, 1.04, 1.03, 1.07, 1.04, 1.03, 1.03, 1.03, 1.04, 1.03,
1.03, 1.04, 1.85, 1.41, 1.35, 1.21, 1.22, 1.15, 1.14, 1.07, 1.1,
1.09, 1.1, 1.09, 1.08, 1.08, 1.09, 1.09, 1.07, 1.06, 1.03, 1.08,
1.05, 1.02, 1.05, 1.99, 1.6, 1.44, 1.4, 1.24, 1.3, 1.21, 1.23,
1.18, 1.18, 1.12, 1.15, 1.09, 1.07, 1.13, 1.1, 1.05, 1.13, 1.09,
1.03, 1.11, 1.07, 1.05), z = structure(c(1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), .Label = c("p1",
"p2", "p3", "p4", "p5"), class = "factor")), class = "data.frame", row.names = c(NA,
-115L))

______________________________________________
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

Fri, Dec 11, 2020 6:25 PM #

On 11/12/2020 6:50 p.m., Rolf Turner wrote:

I don't know of anything easy here.  I think you need to do some tricky 
stuff with formulas to get what you want.  For example,

pred <- function(x, z, ...) {
   a <- unlist(list(...))
   names(a) <- levels(xxx$z)
   x/(x-a[z])
}
strt <- rep(3,length(levels(xxx$z)))
names(strt) <- levels(xxx$z)
fla <- y ~ pred(x, z)
fla[[3]] <- as.call(c(as.list(fla[[3]]), lapply(levels(xxx$z), as.name)))
fit <- nls(fla,start=strt,data=xxx)

That line starting fla[[3]] is the ugly part:  it takes the simple 
formula y ~ pred(x, z) and changes it to y ~ pred(x, z, a1, a2, a3, a4, 
a5).  As far as I know, nls() can't handle vector paramters, it only 
deals with scalars, so this passes them each as a separate argument.

Maybe rlang or one of the other packages like that has code to make this 
less obscure.

Rolf Turner

Fri, Dec 11, 2020 6:39 PM #

On Fri, 11 Dec 2020 19:20:26 -0500

Gabor Grothendieck <ggrothendieck at gmail.com> wrote:

Nya-ha!  I, uh, clearly wasn't thinking clearly!  (Feel a bit *duh*
now!)

Thanks Gabor.

cheers,

Rolf

On Fri, Dec 11, 2020 at 6:51 PM Rolf Turner <r.turner at auckland.ac.nz>
wrote:



I want to fit a model y = x/(x-a) where the value of a depends
on the level of a factor z.  I cannot figure out an appropriate
syntax for nls().  The "parameter" a (to be estimated) should be a
vector of length equal to the number of levels of z.

I tried:

strt <- rep(3,length(levels(z))
names(strt=levels(z)
fit <- nls(y ~ x/(x - a[z]),start=strt,data=xxx)

but of course got an error:

Error in nls(y ~ x/(x - a[z]), start = strt, data = xxx) :
  parameters without starting value in 'data': a

I keep thinking that there is something obvious that I should
be doing, but I can't work out what it is.

Does there *exist* an appropriate syntax for doing what I want
to do?  Can anyone enlighten me?  The data set "xxx" is given
in dput() form at the end of this message.

cheers,

Rolf Turner

--
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

Data set "xxx":

structure(list(x = c(30, 40, 50, 60, 70, 80, 90, 100, 110, 120,
130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250,
30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160,
170, 180, 190, 200, 210, 220, 230, 240, 250, 30, 40, 50, 60,
70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190,
200, 210, 220, 230, 240, 250, 30, 40, 50, 60, 70, 80, 90, 100,
110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230,
240, 250, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140,
150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250), y = c(1.27,
1.16, 1.19, 1.15, 1.09, 1.07, 1.07, 1.05, 1.07, 1.03, 1.05, 1.07,
1.06, 1.03, 1.05, 1.04, 1.03, 1.03, 1.03, 1.02, 1.02, 1.01, 1.01,
1.21, 1.15, 1.1, 1.1, 1.06, 1.06, 1.05, 1.03, 1.07, 1.04, 1.04,
1.02, 1.04, 1.02, 1.04, 1.03, 1.01, 1.03, 1.01, 1, 1.02, 1.03,
1.02, 1.42, 1.27, 1.23, 1.14, 1.17, 1.08, 1.11, 1.06, 1.07, 1.08,
1.06, 1.07, 1.04, 1.03, 1.07, 1.04, 1.03, 1.03, 1.03, 1.04, 1.03,
1.03, 1.04, 1.85, 1.41, 1.35, 1.21, 1.22, 1.15, 1.14, 1.07, 1.1,
1.09, 1.1, 1.09, 1.08, 1.08, 1.09, 1.09, 1.07, 1.06, 1.03, 1.08,
1.05, 1.02, 1.05, 1.99, 1.6, 1.44, 1.4, 1.24, 1.3, 1.21, 1.23,
1.18, 1.18, 1.12, 1.15, 1.09, 1.07, 1.13, 1.1, 1.05, 1.13, 1.09,
1.03, 1.11, 1.07, 1.05), z = structure(c(1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), .Label =
c("p1", "p2", "p3", "p4", "p5"), class = "factor")), class =
"data.frame", row.names = c(NA, -115L))

Andrew Robinson

Fri, Dec 11, 2020 11:39 PM #

Hi Rolf,

It might also be worth experimenting to see if

y ~ 1 / ( 1 - a[z]/x )

yields the same result. It might be cleaner if x appears only once in the expression.

Cheers,

Andrew

--
Andrew Robinson
Director, CEBRA; Professor of Biosecurity
School of BioSciences
University of Melbourne, VIC 3010 Australia
Tel: (+61) 0403 138 955
Email: apro at unimelb.edu.au
Website: http://cebra.unimelb.edu.au/

I acknowledge the traditional owners of the land I inhabit and pay my respects to their elders.

From: R-help <r-help-bounces at r-project.org> on behalf of Rolf Turner <r.turner at auckland.ac.nz>
Sent: Saturday, December 12, 2020 1:39:11 PM
To: Gabor Grothendieck <ggrothendieck at gmail.com>
Cc: r-help at r-project.org <r-help at r-project.org>
Subject: [EXT] Re: [R] nls() syntax

UoM notice: External email. Be cautious of links, attachments, or impersonation attempts


On Fri, 11 Dec 2020 19:20:26 -0500
Gabor Grothendieck <ggrothendieck at gmail.com> wrote:

> The start= argument should be as follows:
>
>  nls(y ~ x/(x - a[z]),start=list(a = strt),data=xxx)

Nya-ha!  I, uh, clearly wasn't thinking clearly!  (Feel a bit *duh*
now!)

Thanks Gabor.

cheers,

Rolf

>
> On Fri, Dec 11, 2020 at 6:51 PM Rolf Turner <r.turner at auckland.ac.nz>
> wrote:
> >
> >
> >
> > I want to fit a model y = x/(x-a) where the value of a depends
> > on the level of a factor z.  I cannot figure out an appropriate
> > syntax for nls().  The "parameter" a (to be estimated) should be a
> > vector of length equal to the number of levels of z.
> >
> > I tried:
> >
> > strt <- rep(3,length(levels(z))
> > names(strt=levels(z)
> > fit <- nls(y ~ x/(x - a[z]),start=strt,data=xxx)
> >
> > but of course got an error:
> >
> > > Error in nls(y ~ x/(x - a[z]), start = strt, data = xxx) :
> > >   parameters without starting value in 'data': a
> >
> > I keep thinking that there is something obvious that I should
> > be doing, but I can't work out what it is.
> >
> > Does there *exist* an appropriate syntax for doing what I want
> > to do?  Can anyone enlighten me?  The data set "xxx" is given
> > in dput() form at the end of this message.
> >
> > cheers,
> >
> > Rolf Turner
> >
> > --
> > Honorary Research Fellow
> > Department of Statistics
> > University of Auckland
> > Phone: +64-9-373-7599 ext. 88276
> >
> > Data set "xxx":
> >
> > structure(list(x = c(30, 40, 50, 60, 70, 80, 90, 100, 110, 120,
> > 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250,
> > 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160,
> > 170, 180, 190, 200, 210, 220, 230, 240, 250, 30, 40, 50, 60,
> > 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190,
> > 200, 210, 220, 230, 240, 250, 30, 40, 50, 60, 70, 80, 90, 100,
> > 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230,
> > 240, 250, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140,
> > 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250), y = c(1.27,
> > 1.16, 1.19, 1.15, 1.09, 1.07, 1.07, 1.05, 1.07, 1.03, 1.05, 1.07,
> > 1.06, 1.03, 1.05, 1.04, 1.03, 1.03, 1.03, 1.02, 1.02, 1.01, 1.01,
> > 1.21, 1.15, 1.1, 1.1, 1.06, 1.06, 1.05, 1.03, 1.07, 1.04, 1.04,
> > 1.02, 1.04, 1.02, 1.04, 1.03, 1.01, 1.03, 1.01, 1, 1.02, 1.03,
> > 1.02, 1.42, 1.27, 1.23, 1.14, 1.17, 1.08, 1.11, 1.06, 1.07, 1.08,
> > 1.06, 1.07, 1.04, 1.03, 1.07, 1.04, 1.03, 1.03, 1.03, 1.04, 1.03,
> > 1.03, 1.04, 1.85, 1.41, 1.35, 1.21, 1.22, 1.15, 1.14, 1.07, 1.1,
> > 1.09, 1.1, 1.09, 1.08, 1.08, 1.09, 1.09, 1.07, 1.06, 1.03, 1.08,
> > 1.05, 1.02, 1.05, 1.99, 1.6, 1.44, 1.4, 1.24, 1.3, 1.21, 1.23,
> > 1.18, 1.18, 1.12, 1.15, 1.09, 1.07, 1.13, 1.1, 1.05, 1.13, 1.09,
> > 1.03, 1.11, 1.07, 1.05), z = structure(c(1L, 1L, 1L, 1L, 1L,
> > 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
> > 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
> > 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
> > 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
> > 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
> > 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
> > 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), .Label =
> > c("p1", "p2", "p3", "p4", "p5"), class = "factor")), class =
> > "data.frame", row.names = c(NA, -115L))

______________________________________________
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.