[R] inverse prediction and Poisson regression

[Spencer Graves]

The problem with nls is that it is NOT maximum likelihood for the 
Poisson distribution.  For the Poisson, the standard deviation is the 
square root of the mean, while nls assumes constant standard deviation. 
  That's why I stayed with glm.  The answers should be comparable, but I 
would prefer the glm answers.  spencer graves

Ravi Varadhan wrote:
> 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
>>
> 
>