From: Ravi Varadhan <rvaradha_at_jhsph.edu>

Date: Thu 14 Jul 2005 - 23:44:47 EST

R-help@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 Received on Thu Jul 14 23:48:11 2005

Date: Thu 14 Jul 2005 - 23:44:47 EST

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@stat.math.ethz.ch [mailto:r-help-
**> bounces@stat.math.ethz.ch] On Behalf Of Stephen B. Cox
**> Sent: Wednesday, July 13, 2005 12:43 PM
**> To: r-help@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@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
*

R-help@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 Received on Thu Jul 14 23:48:11 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:33:40 EST
*