[R] Finney's fiducial confidence intervals of LD50
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 VR'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 EC50Upper
#11.47628 140.8351 8423.78
# while SPSS gives me:
# Lower EC50Upper
#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 EC50Upper
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 EC50Upper
# 46.15749 140.8351 429.7151
# while SPSS gives me:
# Lower EC50Upper
#98,37857 140,83525 205,34483
i.m.s.white [EMAIL PROTECTED] 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
[R] Finney's fiducial confidence intervals of LD50
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: #PRNT50= 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.
Re: [R] Finney's fiducial confidence intervals of LD50
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 D50lowerupper 164.9506 103.3171 191.9777 Best, Renaud 2006/8/21, carlos riveira [EMAIL PROTECTED]: 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: #PRNT50= 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 __ 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.
Re: [R] Finney's fiducial confidence intervals of LD50
carlos riveira carlos.riveira at yahoo.com writes: 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.)? See the example in MASS on dose.p. Don't know if it is comes out the same as SPSS, though. ldose - rep(0:5, 2) numdead - c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16) sex - factor(rep(c(M, F), c(6, 6))) SF - cbind(numdead, numalive = 20 - numdead) budworm.lg0 - glm(SF ~ sex + ldose - 1, family = binomial) dose.p(budworm.lg0, cf = c(1,3), p = 1:3/4) dose.p(update(budworm.lg0, family = binomial(link=probit)), cf = c(1,3), p = 1:3/4) Dieter __ 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.
Re: [R] Finney's fiducial confidence intervals of LD50
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 D50lowerupper 140.8353 103.3171 191.9777 Best, Renaud 2006/8/21, Renaud Lancelot [EMAIL PROTECTED]: 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 D50lowerupper 164.9506 103.3171 191.9777 Best, Renaud 2006/8/21, carlos riveira [EMAIL PROTECTED]: 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: #PRNT50= 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 __ 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.
