[R] Fieller's Conf Limits and EC50's
Ravi Varadhan
rvaradha at jhsph.edu
Thu Jul 14 15:44:47 CEST 2005
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
More information about the R-help
mailing list