[R] inverse prediction and Poisson regression

[Spencer Graves]

	  Peter Dalgaard suggested using "a Michalis-Menten style inhibition 
(mean(Y) = ymax/(1+dose/dose50) or variants thereof)."  If you are not 
familiar with the literature on "Michalis-Menten inhibition", I suggest 
you research that and cite it in your research report.  Even if you 
don't use it, the knowledge of it might help persuade others that you at 
least considered it.

	  If we adopt Ripley's notation, this is mu(dose) = 
ymax/(1+dose/dose50), where ymax and dose50 are to be estimated.  We can 
get this, sort of, using "glm", as follows:  This model is the same as 
log(mu(dose)) = ymax - log(1+dose/dose50).  I have not tried this, but I 
believe that for fixed dose50 = 0.3, this model can be written for glm 
with the toy example I used before as follows:

Phytopath <- data.frame(x=c(0, 0.03, 0.1), 		
	y=c(28, 21, 11))
fitP.3 <- glm(y~offset(-log(1+x/0.3)), data=Phytopath[rep(1:3, 100),],
   family="poisson")
fitP.3$deviance
[1] 359.5893

	  Now, replace 0.3 by other possible values for dose50 to find the 
minimum "deviance".  Then get a 95% confidence interval using "profile 
likelihood" to find the points above and below the minimium where the 
"deviance" exceeds the minimum by  qchisq(0.95, 1) = 3.841459.  There 
are better ways to do this, but none that I can make work in the next 2 
minutes.  When I did this, I get the following:

 > fitP.06 <- glm(y~offset(-log(1+x/0.06)), data=Phytopath[rep(1:3, 100),],
+  family="poisson")
 > fitP.06$deviance
[1] 16.54977
 > fitP.065 <- glm(y~offset(-log(1+x/0.065)), data=Phytopath[rep(1:3, 
100),],
+  family="poisson")
 > fitP.065$deviance
[1] 11.59525
 > fitP.07 <- glm(y~offset(-log(1+x/0.07)), data=Phytopath[rep(1:3, 100),],
+  family="poisson")
 > fitP.07$deviance
[1] 10.96356
 > fitP.075 <- glm(y~offset(-log(1+x/0.075)), data=Phytopath[rep(1:3, 
100),],
+  family="poisson")
 > fitP.075$deviance
[1] 13.49366
 > fitP.08 <- glm(y~offset(-log(1+x/0.08)), data=Phytopath[rep(1:3, 100),],
+  family="poisson")
 > fitP.08$deviance
[1] 18.34463

 From this, I get an approximate maximum likelihood estimate for dose50 
of 0.07 with a 95% confidence interval by rough interpolation of 0.061 
to 0.076.

hope this helps.  spencer graves
p.s.  Have you plotted y vs. x with appropriate transformations?  If 
not, I suggest you do that before anything else.  Doug Bates says that 
students in his classes who don't plot the data get instant F's.  When I 
don't bother plotting the data several different ways, I often do stupid 
things.

Vincent Philion wrote:
> 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,
>