[R] regression for poisson distributed data

[Gabor Grothendieck]

On Tue, Apr 3, 2012 at 11:03 AM, David Winsemius <dwinsemius at comcast.net> wrote:
>
> On Apr 3, 2012, at 9:58 AM, Joachim Audenaert wrote:
>
>> Hello all,
>>
>> I would like to get parameter estimates for different models. For one of
>> them I give the code in example. I am estimating the parameters (i,j and
>> k) with the nls function, which sees the error distribution as normal,
>
>
> Doesn't it really just minimize the squared residuals from a nonlinear
> objective function? I didn't think there was any assumption that there was a
> particular error structure.
>
>
>> I would however like to do the same as nls with the assumption that the
>> errors are poisson distributed.
>
>
>> Is there a way to do this with R? Are there packages designed for this? I
>> tried with the gnm package, but don't understand how to transform my
>> equation to a generalised equation. Is there an option for nls to choose
>> family = poisson?
>
>
>> dat <- read.table(text="N0      FR
>> 10      2
>> 10      3
>> snipped
>>
>> 60      2
>> 60      2
>> 60      5
>> 60      4', header=TRUE)
>
>
> I was wondering if you could just use `glm` in the base/stats package:
>
>
> plot(jitter(FR)~jitter(N0), data=dat)
> ( reg=glm(FR ~ N0, data=dat, family=poisson) )
> lines(10:60, predict(reg, newdata=list(N0=10:60), type="response"))
>
> # It gives you:
>
> FR ~ exp(alpha + beta*N0) + pois-error
>
> # The plot looks adequate. And I would say arguably better than the one I
> get with:
>
>  x <- nls(FR~(exp(i+j*N0)/(1+exp(i+j*N0)))*(k*N0/(k+N0)), data=dat,
>           control =nls.control(maxiter=150, tol=0.01),
>            start=list(i=1,j=.02,k=6))
> summary(x)
> hatx <- predict(x)
> lines(spline(dat$N0,hatx), col="red")
>
> # I also saw Gabor Grothendieck's suggestion which I'm sure is more
> # mathematically sophisticated than my hacking, but when I plot its
> # results, it seems to fit even less well than your original nls function.
>

The normal will tend to overfit to the right side of the plot if the
data is truly poisson since it won't take into account that the
variance of those points is larger.  Note how the normal plot rises
toward the right more than the poisson plot suggesting the normal is
letting those points unduly influence the fit.

-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com