[R] Finney's fiducial confidence intervals of LD50

From: carlos riveira <carlos.riveira_at_yahoo.com>
Date: Tue 22 Aug 2006 - 23:25:02 EST


Thanks for the tip!
  unfortunately though conf limits calculated with Fieller and delta methods do not seem to be in agreement with (and seem to be worse than..) my SPSS results.. Am i doing something wrong?   thanks a lot in advance for your help!!    

  An RSiteSearch on 'Fieller' gave me the following 2 posts   one from S.B.Cox using fieller's conf limits, and one from Varadhan using delta method:    

  1-Fieller's Conf Limits and EC50's
  From: Stephen B. Cox <stephen.cox>
Date: Wed, 13 Jul 2005 11:42:52 -0500        

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

        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) 
 }    

#when i apply it to my data (see below)
    # i get with Fieller's method by Cox

    # Lower EC50 Upper
#11.47628 140.8351 8423.78
# while SPSS gives me:
# Lower EC50 Upper
#98,37857 140,83525 205,34483
   

    library(MASS)
  Response<-c(0,7,26,27,0,5,13,29,0,4,11,25)   Tot<-rep(30.5,12)
  Dose<-rep(c(10,40,160,640),3)
  probit<-glm(formula = Response/Tot~ log10(Dose), family=quasibinomial   (link=probit))
  ec.calc(probit,conf.level=.95,p=.5)
$conf.pts

   Lower EC50 Upper
1.059801 2.148711 3.925507
  $conf.level
[1] 0.95
  $call
ec.calc(obj = probit, conf.level = 0.95, p = 0.5)   Warning message:
multi-argument returns are deprecated in: return(conf.pts, conf.level, call = call)

  2-Fieller's Conf Limits and EC50's
  From: Ravi Varadhan <rvaradha>
Date: Thu, 14 Jul 2005 09:44:47 -0400    

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

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

    #when i apply it to my data
    # i get with delta method by varadhan

    # Lower EC50 Upper
# 46.15749 140.8351 429.7151
# while SPSS gives me:
# Lower EC50 Upper
#98,37857 140,83525 205,34483
  

"i.m.s.white" <i.m.s.white@ed.ac.uk> wrote:   Finney's method for finding the confidence interval for a ratio of parameters is quite simple and is probably described in his book 'Probit analysis'. It is also known as Fieller's method so an RSiteSearch on 'Fieller' might show something useful.

On Mon, Aug 21, 2006 at 08:46:24AM -0700, carlos riveira wrote:

> thanks a lot Renaud. 
> but i was interested in Finney's fiducial confidence intervals of LD50 so to obtain comparable results with SPSS. 
> 
> But your reply leads me to the next question: does anybody know what is the best method (asymptotic, bootstrap etc.) for calculating confidence intervals of LD50? 
> 
> i could "get rid" of Finney's fiducial confidence intervals but only if there was a better method..
> 
> any idea? 
> 
> 
> 
> Renaud Lancelot wrote:
> Date: Mon, 21 Aug 2006 16:35:49 +0200
> From: "Renaud Lancelot" 
> To: "carlos riveira" 
> Subject: Re: [R] Finney's fiducial confidence intervals of LD50
> CC: r-help@stat.math.ethz.ch
> 
> Sorry there was a typo in my previous reply:
> 
> > D50 <- 10^c(logD50 + c(0, -1.96, 1.96) * attr(logD50, "SE"))
> > names(D50) <- c("D50", "lower", "upper")
> > D50
> D50 lower upper
> 140.8353 103.3171 191.9777
> 
> Best,
> 
> Renaud
> 
> 2006/8/21, Renaud Lancelot :
> > I don't know what Finney's fiducial confidence interval is but if your
> > problem is to handle the output of dose.p (from MASS), you can do as
> > follows:
> >
> > > library(MASS)
> > > Response <- c(0, 7, 26, 27, 0, 5, 13, 29, 0, 4, 11, 25)
> > > Tot <- rep(30.5, 12)
> > > Dose <- rep(c(10, 40, 160, 640), 3)
> > > fm <- glm(Response/Tot ~ log10(Dose), family = quasibinomial(link = probit))
> > > logD50 <- dose.p(fm, cf = 1:2, p = 0.5)
> > > D50 <- 10^c(logD50 + c(1, -1.96, 1.96) * attr(logD50, "SE"))
> > > names(D50) <- c("D50", "lower", "upper")
> > > D50
> > D50 lower upper
> > 164.9506 103.3171 191.9777
> >
> > Best,
> >
> > Renaud
> >
> > 2006/8/21, carlos riveira :
> > > I am working with Probit regression (I cannot switch to logit) can anybody help me in finding out how to obtain with R Finney's fiducial confidence intervals for the levels of the predictor (Dose) needed to produce a proportion of 50% of responses(LD50, ED50 etc.)?
> > > If the Pearson chi-square goodness-of-fit test is significant (by default), a heterogeneity factor should be used to calculate the limits.
> > >
> > > Response<-c(0,7,26,27,0,5,13,29,0,4,11,25)
> > > Tot<-rep(30.5,12)
> > > Dose<-rep(c(10,40,160,640),3)
> > > probit<-glm(formula = Response/Tot~ log10(Dose), family=quasibinomial
> > > (link=probit))
> > > D50<- round(10^(dose.p(probit,cf=1:2,p=0.5)))
> > >
> > > #This is what SPSS calculates. I would like to reproduce these results with R:
> > > #SPSS RESULTS:
> > > #D50= 140,83525
> > > #CI = 98,37857;205,34483
> > > #Regr.coeff= 1,91676 (S.E.=0,16765)
> > > #Intercept=-4,11856 (S.E.=0,36355)
> > > Thanks a lot for your help.
> > >
> > > Carlos
> > >
> > > __________________________________________________
> > >
> > >
> > >
> > > [[alternative HTML version deleted]]
> > >
> > > ______________________________________________
> > > 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
> > > and provide commented, minimal, self-contained, reproducible code.
> > >
> >
> >
> > --
> > Renaud LANCELOT
> > D?partement Elevage et M?decine V?t?rinaire (EMVT) du CIRAD
> > Directeur adjoint charg? des affaires scientifiques
> >
> > CIRAD, Animal Production and Veterinary Medicine Department
> > Deputy director for scientific affairs
> >
> > Campus international de Baillarguet
> > TA 30 / B (B?t. B, Bur. 214)
> > 34398 Montpellier Cedex 5 - France
> > T?l +33 (0)4 67 59 37 17
> > Secr. +33 (0)4 67 59 39 04
> > Fax +33 (0)4 67 59 37 95
> >
> 
> 
> -- 
> Renaud LANCELOT
> D?partement Elevage et M?decine V?t?rinaire (EMVT) du CIRAD
> Directeur adjoint charg? des affaires scientifiques
> 
> CIRAD, Animal Production and Veterinary Medicine Department
> Deputy director for scientific affairs
> 
> Campus international de Baillarguet
> TA 30 / B (B?t. B, Bur. 214)
> 34398 Montpellier Cedex 5 - France
> T?l +33 (0)4 67 59 37 17
> Secr. +33 (0)4 67 59 39 04
> Fax +33 (0)4 67 59 37 95
> 
> 
> 
> ---------------------------------
> 
> [[alternative HTML version deleted]]
> 

> ______________________________________________
> 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 > and provide commented, minimal, self-contained, reproducible code.
-- 
************************************************
* I.White *
* University of Edinburgh *
* Ashworth Laboratories, West Mains Road *
* Edinburgh EH9 3JT *
* Fax: 0131 650 6564 Tel: 0131 650 5490 *
* E-mail: i.m.s.white@ed.ac.uk *
************************************************














 		
---------------------------------
Stay in the know. Pulse on the new Yahoo.com.  Check it out. 
	[[alternative HTML version deleted]]

______________________________________________
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
and provide commented, minimal, self-contained, reproducible code.
Received on Wed Aug 23 00:52:23 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Wed 23 Aug 2006 - 02:22:14 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.