Hi,
In a theoretical experiment, I expose individual subjects to various doses of a toxic agent (example here is radiation), and later determine if the subject is dead or alive.
I want to use either logit or probit to describe probability of death as a function of dose.
Specifically, I want to calculate the dose that will kill 50% of the subjects (LD50), and I want to report LD50 with 95% CI.
Finally, I want to test if the logit LD50 is different from the probit LD50.
The last question doesn't make much sense, but if I understand how to do the test, I can use it on more interesting questions (like, if I only use logit and give a radio protector to some subjects, will they have a significantly higher LD50?).
My questions are:
a) How do I calculate SE on LD50, and use that to calculate 95% CI? I have seen this question raised in various forums and publications. It seems complicated, and I don't know which formula to use, or how to enter it in Stata in a way that refers only to stored results (like _b[cons]).
b) What is the best reference for how to calculate SE?
c) Have I made any other mistakes in the code below, including how to test the difference in LD50s and for the plot used for visual inspection of the data and the logit/probit functions?
d) This is my first post, and please correct any mistakes I may have done on how to post a question on this list.
Any help is greatly appreciated,
Many thanks,
Jan
In a theoretical experiment, I expose individual subjects to various doses of a toxic agent (example here is radiation), and later determine if the subject is dead or alive.
I want to use either logit or probit to describe probability of death as a function of dose.
Specifically, I want to calculate the dose that will kill 50% of the subjects (LD50), and I want to report LD50 with 95% CI.
Finally, I want to test if the logit LD50 is different from the probit LD50.
The last question doesn't make much sense, but if I understand how to do the test, I can use it on more interesting questions (like, if I only use logit and give a radio protector to some subjects, will they have a significantly higher LD50?).
My questions are:
a) How do I calculate SE on LD50, and use that to calculate 95% CI? I have seen this question raised in various forums and publications. It seems complicated, and I don't know which formula to use, or how to enter it in Stata in a way that refers only to stored results (like _b[cons]).
b) What is the best reference for how to calculate SE?
c) Have I made any other mistakes in the code below, including how to test the difference in LD50s and for the plot used for visual inspection of the data and the logit/probit functions?
d) This is my first post, and please correct any mistakes I may have done on how to post a question on this list.
Code:
. version 14
. webuse auto, clear // Generate data test set
(1978 Automobile Data)
. generate subjectID=_n // Clarify that observations are individual subjects
. generate dose=(7500-weight)/100 // Exposure (radiation dose)
. label variable dose "Radiation dose (Gy)"
. generate survival=foreign // Outcome (survival, 0 or 1)
. label variable survival "Survival: 0=alive 1=dead"
. drop w f m* p r h t* l di g
. logit survival dose // logit
Iteration 0: log likelihood = -45.03321
Iteration 1: log likelihood = -30.669507
Iteration 2: log likelihood = -29.068209
Iteration 3: log likelihood = -29.054006
Iteration 4: log likelihood = -29.054002
Iteration 5: log likelihood = -29.054002
Logistic regression Number of obs = 74
LR chi2(1) = 31.96
Prob > chi2 = 0.0000
Log likelihood = -29.054002 Pseudo R2 = 0.3548
------------------------------------------------------------------------------
survival | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
dose | .2587387 .0609368 4.25 0.000 .1393048 .3781727
_cons | -13.12281 3.018091 -4.35 0.000 -19.03816 -7.207457
------------------------------------------------------------------------------
. scalar _consl = _b[_cons] // Used for plot
. scalar dosel = _b[dose] // Used for plot
. scalar LD50l = -_b[_cons]/_b[dose] // LD50 (logit)
. *scalar LD50lse = ??? // SE - Don't know how to do this
. scalar LD50lse = 5 // Just a random number to run the next part
. scalar LD50lupp = LD50l+(1.96*LD50lse) // Don't know if can be done on this scale
. scalar LD50llow = LD50l-(1.96*LD50lse)
. display "LD50(logit): " %4.2f LD50l " Gy"
LD50(logit): 50.72 Gy
. display "95% CI: " %4.2f LD50llow "-" %4.2f LD50lupp
95% CI: 40.92-60.52
. probit survival dose // probit
Iteration 0: log likelihood = -45.03321
Iteration 1: log likelihood = -29.534424
Iteration 2: log likelihood = -28.912832
Iteration 3: log likelihood = -28.908407
Iteration 4: log likelihood = -28.908407
Probit regression Number of obs = 74
LR chi2(1) = 32.25
Prob > chi2 = 0.0000
Log likelihood = -28.908407 Pseudo R2 = 0.3581
------------------------------------------------------------------------------
survival | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
dose | .15049 .0326453 4.61 0.000 .0865064 .2144735
_cons | -7.631124 1.602292 -4.76 0.000 -10.77156 -4.490688
------------------------------------------------------------------------------
. scalar _consp = _b[_cons] // Used for plot
. scalar dosep = _b[dose] // Used for plot
. scalar LD50p = -_b[_cons]/_b[dose] // LD50 (probit)
. *scalar LD50pse = ??? // SE - Don't know how to do this
. scalar LD50pse = 5 // Just a random number to run the next part
. scalar LD50pupp = LD50p+(1.96*LD50pse) // Don't know if can be done on this scale
. scalar LD50plow = LD50p-(1.96*LD50pse)
. display "LD50(probit): " %4.2f LD50p " Gy"
LD50(probit): 50.71 Gy
. display "95% CI: " %4.2f LD50plow "-" %4.2f LD50pupp
95% CI: 40.91-60.51
. * Test difference between LD50logit and LD50probit
. scalar Z = (LD50l-LD50p)/sqrt(LD50lse^2+LD50pse^2)
. scalar Pval = 2*normal(-abs(Z))
. display "P value:" %9.2e Pval
P value: 9.99e-01
. * Plot data and logit/probit functions
. twoway function y = invlogit(_consl + dosel*x), range(dose) || ///
> function y = normprob(_consp + dosep*x), range(dose) || ///
> scatter survival dose, ///
> xtitle("Radiation dose (Gy)") ytitle("Probability of death") ///
> legend(order(1 2) label(1 "logit") label(2 "probit"))
.
Any help is greatly appreciated,
Many thanks,
Jan
Comment