[R] regression for poisson distributed data

[David Winsemius]

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.

attach(dat)  #  I don't normally use attach but it seemed quicker at  
that moment.
Mean <- function(p) { i <- p[1]; j <- p[2]; k <- p[3];
exp(i+j*N0)/(1+exp(i+j*N0))*(k*N0/(k+N0)) }
ll <- function(p) - sum(dpois(FR, Mean(p), log = TRUE))
out <- optim(c(1, 1, 1), ll); out

lines(Mean(out$par) ~ N0, col="blue")
  
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Rplot.pdf
Type: application/pdf
Size: 84147 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20120403/082e0133/attachment.pdf>
-------------- next part --------------


All my comments are moot if there is strong theory that demands this  
functional form with Poisson errors, but you have not yet provided  
background that suggest this to be the case.


>
> Lower in the mail the code with the model and visualisations I use to
> check my results.

Nope. Attachments are scrubbed if not an acceptable type to the  
server. (PDF's are accepted.)

> I also copied the test dataset from my txt file.
>
> I am using R 2.15 and Rstudio to visualise it.
>
>
> plot(FR~N0)
> x <-
> nls(FR~(exp(i+j*N0)/(1+exp(i+j*N0)))*(k*N0/(k 
> +N0)),start=list(i=0.02,j=0.002,k=6))
> summary(x)
> hatx <- predict(x)
> lines(spline(N0,hatx))
>
> Kind Regards
> Joachim
>
> Don't waste paper! Think about the environment...
> ____________________________________
> 	[[alternative HTML version deleted]]

Don't waste electrons ( or our time in accessing the Archives). Think  
about our mailing list environment and the Posting Guide which tells  
you NOT to use HTML.

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

David Winsemius, MD
West Hartford, CT