[R] how to code y~x/(x+a) in lm() function

Rolf Turner rolf.turner at xtra.co.nz
Fri Aug 23 01:12:16 CEST 2013


I tried your suggestion on the toy example that was posted by David Carlson
a while back:

set.seed(42)
x <- 1:15
y <- x/(1+x) + rnorm(15,0,0.02)

I found that I had to:

     (a) Wrap the "1/x" inside "I()".
     (b) Set the offset term to be rep(1,length(x)).

Bottom line:

fit <- glm(y ~ 
0+I(1/x)+offset(rep(1,length(x))),family=gaussian(link="inverse")
plot(x,y)
lines(x,fitted(fit))

seems to give a "sensible" picture.

Doing

nlsfit <- nls(y ~ x/(a+x),start=list(a=0.9))
lines(x,fitted(nlsfit),col="red")

indicates that the two fits are "visually indistinguishable". Both 
methods estimate
"a" as 1.061.

     cheers,

         Rolf

On 22/08/13 12:54, Ben Bolker wrote:
> 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?



More information about the R-help mailing list