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: [EMAIL PROTECTED] [mailto:r-help- > [EMAIL PROTECTED] 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