[R] inverse prediction and Poisson regression

[Peter Dalgaard BSA]

Prof Brian Ripley <ripley at stats.ox.ac.uk> writes:


> > x	y
> > 0		28
> > 0.03		21
> > 0.1		11
> > 0.3		15
> > 1		5
> > 3		4
> > 10		1
> > 30		0
> > 100		0
....
> > Where X is dose and Y is response. 
> > the relation is linear for log(response) = b log(dose) + intercept
> 
> Is that log(*mean* response), that is a log link and exponential decay 
> with dose?
> 
> > Response for dose 0 is a "control" = Ymax. So, What I want is the dose
> > for 50% response. For instance, in example 1:
> > 
> > Ymax = 28 (this is also an observation with Poisson error)
> 
> Once you observe Ymax, Y is no longer Poisson.
> 
> > So I want dose for response = 14 = approx. 0.3
> 
> What exactly is Ymax?  Is it the response at dose 0?  The mean response at
> dose 0?  The largest response?  About the only thing I can actually
> interpret is that you want to fit a curve of mean response vs dose, and
> find the dose at which the mean response is half of that at dose 0.
> That one is easy.
> 
> I think you are confusing response with mean response, and we can't 
> disentangle them for you.

I don't feel all that confused. Y is Poisson distributed with some
mean depending on x. Ymax is a value at X=0, i.e. Poisson distr.
with a mean as large as it can be. 

I think the main confusion here is trying to fit a functional
relationship which doesn't extend to X=0. If you extrapolate a
log-loglinear relation back to X=0, you get an infinite maximal
response if b is negative, so this is going to be inconsistent with a
finite Ymax. In some of the data sets I believe you actually do see a
leveling off for very small doses.

If you insist on this peculiar model, you'd end up with estimating the
mean of Ymax by its observed value. Then you can get b and the
intercept from the observations with X>0, and find your estimate of
halving dose by solving

 log(Ymax/2) = b * log(dose50) + intercept

i.e. dose50 = (log(Ymax/2)-intercept)/b. That's a nonlinear function
of the estimates, so you'd need (at least) to employ the Delta method
to find the approximate variance of the estimate.

However, I'd suggest that you should look for a more realistic
functional form of the relation, e.g. a logistic curve in log(x) or a
Michalis-Menten style inhibition (mean(Y) = ymax/(1+dose/dose50) or
variants thereof). These models are not (necessarily) GLMs, but I
think you can fit them quite well with gnls() and a suitable variance
function.

[In fact the ymax/(1+dose/dose50) model is a GLM if you use an inverse
link and reparametrize with 1/ymax, 1/(ymax*dose50), but inverse links
are not built-in for the poisson() family, so you'd have to modify the
code yourself.]

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907