[R] inverse prediction and Poisson regression

[Ravi Varadhan]

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
>