Hey All, I wanna to fit a model y~x/(a+x) to my data, here is the code I use now: lm((1/y-1)~I(1/x)+0, data=b) and it will return the coefficient which is value of a however, if I use the code above, I am not able to draw a curve the presents this equation. How can I do this? Thanks for your help! [[alternative HTML version deleted]]
(1) It is not acceptable to use "wanna" in written English. You
should say
"I want to fit a model ....".
(2) The model you have fitted is *not* equivalent to the model you first
state.
If you write "y ~ x/(a+x)" you are tacitly implying that
y = x/(a+x) + E
where the "errors" E are i.i.d. with mean 0.
If this is the case then it will *not* be the case that
1/y = 1 + a/x + E
with the E values being i.i.d. with mean 0.
If the model "y ~ x/(a+x)" is really what you want to fit, then you
should
be using non-linear methods, e.g. by applying the function nls().
cheers,
Rolf Turner
On 21/08/13 09:39, Ye Lin wrote:> Hey All,
>
> I wanna to fit a model y~x/(a+x) to my data, here is the code I use now:
>
> lm((1/y-1)~I(1/x)+0, data=b)
>
> and it will return the coefficient which is value of a
>
> however, if I use the code above, I am not able to draw a curve the
> presents this equation. How can I do this?
?curve set.seed(42) x <- 1:15 y <- x/(1+x)+rnorm(15, 0, .02) plot(y~x) lm.out <- lm((1/y-1)~I(1/x)+0) curve(x/(coef(lm.out)+x), 1, 15, add=TRUE) ------------------------------------- David L Carlson Associate Professor of Anthropology Texas A&M University College Station, TX 77840-4352 -----Original Message----- From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Ye Lin Sent: Tuesday, August 20, 2013 4:40 PM To: R help Subject: [R] how to code y~x/(x+a) in lm() function Hey All, I wanna to fit a model y~x/(a+x) to my data, here is the code I use now: lm((1/y-1)~I(1/x)+0, data=b) and it will return the coefficient which is value of a however, if I use the code above, I am not able to draw a curve the presents this equation. How can I do this? Thanks for your help! [[alternative HTML version deleted]] ______________________________________________ R-help at r-project.org mailing list 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.
On 13-08-21 05:17 PM, Rolf Turner wrote:> > > Thott about this a bit more and have now decided that I don't understand > after all. > > Doesn't > > glm(1/y~x,family=gaussian(link="inverse")) > > fit the model > > 1/y ~ N(1/(a+bx), sigma^2) > > whereas what the OP wanted was > > y ~ N(x/(a+x),sigma^2) ???I goofed slightly. y ~ 1/x with inverse link gives 1/y = a + b*(1/x) y = 1/(a+b*(1/x)) = x/(a*x+b) Hmmm. Is there an offset trick we can use? y = x/(a+x) 1/y = (a+x)/x 1/y = (a/x) + 1 1/y = a*(1/x) + 1 So I *think* we want glm(y~1/x+0+offset(1),family=gaussian(link="inverse")) I'm forwarding this back to r-help.> > I can't see how these models can be equivalent. What am I missing? > > cheers, > > Rolf > > > > On 22/08/13 03:49, Ben Bolker wrote: >> Rolf Turner <rolf.turner <at> xtra.co.nz> writes: >> >>> On 21/08/13 11:23, Ye Lin wrote: >>>> T >>>> hanks for your insights Rolf! The model I want to fit is y=x/a+x with >>>> no intercept, so I transformed it to 1/y=1+a/x as they are the same. >>> For crying out loud, they are ***NOT*** the same. The equations y >>> x/(a+x) and >>> 1/y = 1 + a/x are indeed algebraically identical, but if an "error" or >>> "noise" term is added >>> to each then then the nature of the error term is vastly different. It >>> is the error or >>> noise term that is of central concern in a statistical context. >>> >>> cheers, >>> >>> Rolf >> >> For what it's worth this model can also be fitted (without messing >> up the error structure) via >> >> glm(1/y~x,family=gaussian(link="inverse")) >> >> Although you may not get the parameters in exactly the form you >> want. >> >> ______________________________________________ >> R-help at r-project.org mailing list >> 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. >> >