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

[Ravi Varadhan]

Hi,

I didn't verify your formulas for Fieller's method of computing the
confidence interval. A slightly simpler approach is to use the Delta method
to compute the CI.  It is also valid for any link function.  It yields a
simpler formula for the variance of EC50 (for any link function):

varEC50 <- 1/b1^2 * (var.b0 + EC50^2*var.b1 + 2*EC50*cov.b0.b1)

So, you can compute the CIs as:

LCL <- EC50 - zalpha.2 * sqrt(varEC50)
UCL <- EC50 + zalpha.2 * sqrt(varEC50)

This works for any EC_p, where p is the probability of getting a positive
response, and for link function.  The CI from the Delta method should be
very nearly the same as that obtained using Fieller's method for EC50.  For
smaller probabilities (e.g., p < 0.1), CIs for the EC_p values obtained
using the two methods can be slightly different.

Ravi.

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-
> bounces at stat.math.ethz.ch] On Behalf Of Stephen B. Cox
> Sent: Wednesday, July 13, 2005 12:43 PM
> To: r-help at stat.math.ethz.ch
> Subject: [R] Fieller's Conf Limits and EC50's
> 
> 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
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-
> guide.html