[R] Fieller's Conf Limits and EC50's

[Stephen B. Cox]

Folks

I have modified an existing function to calculate 'ec/ld/lc' 50 values 
and their associated Fieller's confidence limits.  It is based on 
EC50.calc (writtien by John Bailer)  - but also borrows from the dose.p 
(MASS) function.  My goal was to make the original EC50.calc function 
flexible with respect to 1) probability at which to calculate the 
expected dose, and 2) the link function.  I would appreciate comments 
about the validity of doing so!  In particular - I want to make sure 
that the confidence limit calculations are still valid when changing the 
link function.

ec.calc<-function(obj,conf.level=.95,p=.5) {

 # calculates confidence interval based upon Fieller's thm.
 # modified version of EC50.calc found in P&B Fig 7.22
 # now allows other link functions, using the calculations
 # found in dose.p (MASS)
 # SBC 19 May 05

        call <- match.call()

         coef = coef(obj)
         vcov = summary.glm(obj)$cov.unscaled
         b0<-coef[1]
         b1<-coef[2]
         var.b0<-vcov[1,1]
         var.b1<-vcov[2,2]
         cov.b0.b1<-vcov[1,2]
         alpha<-1-conf.level
         zalpha.2 <- -qnorm(alpha/2)
         gamma <- zalpha.2^2 * var.b1 / (b1^2)
         eta = family(obj)$linkfun(p)  #based on calcs in V&R's dose.p

         EC50 <- (eta-b0)/b1

         const1 <- (gamma/(1-gamma))*(EC50 + cov.b0.b1/var.b1)

         const2a <- var.b0 + 2*cov.b0.b1*EC50 + var.b1*EC50^2 -
                    gamma*(var.b0 - cov.b0.b1^2/var.b1)

         const2 <- zalpha.2/( (1-gamma)*abs(b1) )*sqrt(const2a)

         LCL <- EC50 + const1 - const2
         UCL <- EC50 + const1 + const2

         conf.pts <- c(LCL,EC50,UCL)
         names(conf.pts) <- c("Lower","EC50","Upper")

         return(conf.pts,conf.level,call=call)
 }


Thanks

Stephen Cox