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, 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? Lower in the mail the code with the model and visualisations I use to check my results. 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)) N0 FR 10 2 10 3 10 2 10 4 10 2 10 2 10 1 10 2 10 2 10 2 20 2 20 3 20 3 20 3 20 4 20 2 20 4 20 2 20 3 20 2 30 1 30 2 30 3 30 4 30 5 30 6 30 2 30 3 30 2 30 2 40 2 40 3 40 3 40 6 40 5 40 4 40 3 40 3 40 2 40 3 50 4 50 5 50 2 50 3 50 7 50 5 50 4 50 3 50 4 50 5 60 5 60 6 60 8 60 4 60 4 60 3 60 2 60 2 60 5 60 4 Kind Regards Joachim Don't waste paper! Think about the environment before printing this e-mail ______________________________________ Joachim Audenaert Adviesdienst Gewasbescherming Proefcentrum voor Sierteelt Schaessestraat 18 B-9070 Destelbergen Tel. +32 9 353 94 71 Fax +32 9 353 94 95 E-mail: joachim.audenaert@pcsierteelt.be www.pcsierteelt.be ______________________________________ [[alternative HTML version deleted]]
On Tue, Apr 3, 2012 at 9:58 AM, Joachim Audenaert <Joachim.Audenaert at pcsierteelt.be> 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, 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? > > Lower in the mail the code with the model and visualisations I use to > check my results. 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)) >Just minimize the negative log likelihood from first principles. Define the mean function and then define the negative log likelihood in terms of that: 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 plot(FR ~ N0, pch = 20) lines(Mean(out$par) ~ N0) -- Statistics & Software Consulting GKX Group, GKX Associates Inc. tel: 1-877-GKX-GROUP email: ggrothendieck at gmail.com
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
?Just write down the loglikelihood function and send it to optim? Kjetil On Tue, Apr 3, 2012 at 8:58 AM, Joachim Audenaert <Joachim.Audenaert at pcsierteelt.be> 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, 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? > > Lower in the mail the code with the model and visualisations I use to > check my results. 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)) > > N0 ? ? ?FR > 10 ? ? ?2 > 10 ? ? ?3 > 10 ? ? ?2 > 10 ? ? ?4 > 10 ? ? ?2 > 10 ? ? ?2 > 10 ? ? ?1 > 10 ? ? ?2 > 10 ? ? ?2 > 10 ? ? ?2 > 20 ? ? ?2 > 20 ? ? ?3 > 20 ? ? ?3 > 20 ? ? ?3 > 20 ? ? ?4 > 20 ? ? ?2 > 20 ? ? ?4 > 20 ? ? ?2 > 20 ? ? ?3 > 20 ? ? ?2 > 30 ? ? ?1 > 30 ? ? ?2 > 30 ? ? ?3 > 30 ? ? ?4 > 30 ? ? ?5 > 30 ? ? ?6 > 30 ? ? ?2 > 30 ? ? ?3 > 30 ? ? ?2 > 30 ? ? ?2 > 40 ? ? ?2 > 40 ? ? ?3 > 40 ? ? ?3 > 40 ? ? ?6 > 40 ? ? ?5 > 40 ? ? ?4 > 40 ? ? ?3 > 40 ? ? ?3 > 40 ? ? ?2 > 40 ? ? ?3 > 50 ? ? ?4 > 50 ? ? ?5 > 50 ? ? ?2 > 50 ? ? ?3 > 50 ? ? ?7 > 50 ? ? ?5 > 50 ? ? ?4 > 50 ? ? ?3 > 50 ? ? ?4 > 50 ? ? ?5 > 60 ? ? ?5 > 60 ? ? ?6 > 60 ? ? ?8 > 60 ? ? ?4 > 60 ? ? ?4 > 60 ? ? ?3 > 60 ? ? ?2 > 60 ? ? ?2 > 60 ? ? ?5 > 60 ? ? ?4 > > Kind Regards > Joachim > > Don't waste paper! Think about the environment before printing this e-mail > > ______________________________________ > > Joachim Audenaert > Adviesdienst Gewasbescherming > Proefcentrum voor Sierteelt > Schaessestraat 18 > B-9070 Destelbergen > Tel. +32 9 353 94 71 > Fax +32 9 353 94 95 > E-mail: joachim.audenaert at pcsierteelt.be > www.pcsierteelt.be > ______________________________________ > ? ? ? ?[[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org mailing list > 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.