[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