regression for poisson distributed data
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