[R] inverse prediction and Poisson regression
Ravi Varadhan
rvaradha at jhsph.edu
Fri Jul 25 17:55:22 CEST 2003
Vincent:
Here is a simple solution using Prof. Bates' non-linear least squares
algorithm:
Best,
Ravi.
> Phytopath <- data.frame(x=c(0, 0.03, 0.1), y=c(28, 21, 11))
> Phyto.nls <- nls(y ~ Ymax/(1 + x/x50),data=Phytopath,start=list
(Ymax=20.0,x50=0.01),trace=T)
404.3058 : 20.00 0.01
15.76932 : 27.96313636 0.04960484
2.043625 : 28.2145584 0.0694645
1.851401 : 28.33886844 0.07198951
1.851231 : 28.34892493 0.07185953
1.851230 : 28.34843670 0.07186804
1.851230 : 28.3484688 0.0718675
> summary(Phyto.nls)
Formula: y ~ Ymax/(1 + x/x50)
Parameters:
Estimate Std. Error t value Pr(>|t|)
Ymax 28.34847 1.31522 21.554 0.0295 *
x50 0.07187 0.01348 5.332 0.1180
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 1.361 on 1 degrees of freedom
Correlation of Parameter Estimates:
Ymax
x50 -0.6001
----- Original Message -----
From: Vincent Philion <vincent.philion at irda.qc.ca>
Date: Friday, July 25, 2003 9:25 am
Subject: Re: [R] inverse prediction and Poisson regression
> Hi, ... and good morning!
>
> ;-)
>
> On 2003-07-25 08:43:35 -0400 Spencer Graves
> <spencer.graves at PDF.COM> wrote:
>
> > The Poisson assumption means that Y is a number of
> independent events from
> > a theoretically infinite population occurring in a specific time
> or place.
> > The function "glm" with 'family="poisson"' with the default link
> = "log"
> > assumes that the logarithm of the mean of Y is a linear model in
> the
> > explanatory variable.
>
> OK, I think my data can fit that description.
>
> >
> > How is Y measured?
>
> Y is the number of line intercepts which encounters mycelial
> growth. i/e if mycelia intercepts the line twice, 2 is reported.
> This follows poisson.
>
> If it the number out N, with N approximately 500 (and you know N),
> > then you have a logistic regression situation.
>
> No, 500 spores can grow, but there is no "real" limit on the
> amount of growth possible, and so no limit on the number of
> intercepts. So this is why I adopted Poisson, not knowing how
> complicated my life would become!!!
> ;-)
>
> In that case, section 7.2 in
> > Venables and Ripley (2002) should do what you want. If Y is a
> percentage
> > increase
>
> ... But you may be right, that I'm making this just too
> complicated and that I should simply look at percentage... Any
> comments on that?
>
>
> > When dose = 0, log(dose) = (-Inf). Since 0 is a legitimate
> dose,
> > log(dose) is not acceptable in a model like this. You need a
> model like
> > Peter suggested.
>
> OK, I see I will need stronger coffee to tackle this, but I will
> read this in depth today.
>
> Depending on you purpose, log(dose+0.015) might be
> > sufficiently close to a model like what Peter suggested to
> answer your
> > question. If not, perhaps this solution will help you find a
> better
> > solution.
>
> In other words, "cheat" and model Y_0 with a "small" value =
> log(0.015) ? How would this affect the LD50 value calculated and
> the confidence intervals? I guess I could try several methods, but
> how would I go about choosing the right one? Criteria?
>
> > I previously was able to get dose.p to work in R, and I just
> now was able
> > to compute from its output. The following worked in both S-Plus
> 6.1 and R
> > 1.7.1:
> >
> >> LD50P100p <- print(LD50P100)
> > Dose SE
> > p = 14: -2.451018 0.04858572
> >> exp(LD50P100p[1,1]+c(-2,0,2)*LD50P100p[1,2])-0.015
> > [1] 0.06322317 0.07120579 0.08000303
>
> OK, I will need to try this (later today). I don't see "dose.p" in
> this?
> again, many thanks,
>
> --
> Vincent Philion, M.Sc. agr.
> Phytopathologiste
> Institut de Recherche et de Développement en Agroenvironnement (IRDA)
> 3300 Sicotte, St-Hyacinthe
> Québec
> J2S 7B8
>
> téléphone: 450-778-6522 poste 233
> courriel: vincent.philion at irda.qc.ca
> Site internet : www.irda.qc.ca
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>
More information about the R-help
mailing list