[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

From |
"Rosy Reynolds" <[email protected]> |

To |
"statalist" <[email protected]> |

Subject |
Re: st: reverse prediction - confidence interval for x at given y in nonlinear model |

Date |
Thu, 25 Oct 2007 14:37:22 +0100 |

Dear Maarten,

That's absolutely brilliant, thank you very much. The general formula was

the thing I had failed to find, so that is extremely useful, and it is

always helpful to see example code too. For any arbitrary y, I can now

recast it as an ED_{100*a}, and work from there.

Thanks also for the warning about extrapolating beyond our observed x. That

is not such a big problem for us, in practice, because we are working with

bacteria in vitro, so we don't have the same tolerability/toxicity/safety

issues to consider as we would if we were dosing people or animals (or even

plants). We can use very wide dose ranges indeed.

Apologies for the delayed reply. My e-mail system is erratic about sending today.

best wishes

Rosy

BSAC Resistance Surveillance Co-ordinator

www.bsacsurv.org

P.S. Thanks too, in general, to the many regular contributors to Statalist.

I have learned a lot by reading your responses to other people's questions,

and I appreciate so many of you contributing your time and expertise as you

do.

----- Original Message ----- From: "Maarten buis" <[email protected]>

To: <[email protected]>

Sent: Thursday, October 25, 2007 12:30 PM

Subject: Re: st: reverse prediction - confidence interval for x at given y

in nonlinear model

--- Rosy Reynolds <[email protected]> wrote:Sigmoid models are customary in pharmacodynamics (dose-response studies). According to custom, I am using a 4-parameter logistic (sigmoid Emax) model. This is very easily done with -nl- as Stata has this model already built in. The model is y= b0 + b1/(1 + exp(-b2*(x-b3))) + error and the coefficients can be interpreted as b0 = baseline outcome b1 = Emax i.e. largest change from baseline b2 = Hill or slope coefficient b3 = ED50 i.e. value of x (dose) required to produce half-maximal effect, that is x required for y=b0 + b1 / 2 As the ED50 is a parameter of the model, -nl- reports it with a standard error and confidence interval. What I would like to do is obtain estimates with standard errors and confidence intervals for other similar measures e.g. the ED90, the dose required for 90% of maximal effect. In general, how could I calculate an estimate and confidence interval for the x required to achieve any given value of y?The general formula for a ED_{100*a} = ln(a/(1-a))/b2 + b3. You can use -nlcom- to calculate ED for any a, with its confidence interval. However, you should be aware that the estimated ED's can be way outside your observed range of x, in other words you can end up with an extrapolation. In the example below I estimated different EDs and also, just for fun, created a graph of ED_{100*a} against a. In this graph I added a rug for the values of x, and two horizontal lines representing the minimum and maximum observed value of x, to show when these estimated EDs are based on sparse data and when they are actually extrapolations. I created the rug by putting a minor tick on the inside of the y-axis for each unique value of x. The unique values of x are stored in r(levels) after -levelsof(x)-, and the minor ticks were placed using the -ymticks()- option. Tricks I used to store the lower and upper bounds of the confidence interval are discussed in: http://home.fsw.vu.nl/m.buis/wp/pvalue.html *---------------------- begin example -------------------- /*Create some data*/ set more off set seed 1234 drop _all set obs 500 gen x = invnorm(uniform()) gen y = 2 + 4 /(1 + exp(-.5*(x - 1))) + .5*invnorm(uniform()) /*Estimate the model*/ nl log4: y x estimates store log4 /*Calculating ED90*/ local lodds = ln(.9/(1-.9)) nlcom `lodds'/[b2]_b[_cons] + [b3]_b[_cons] /*graphing EDs*/ gen a = .05*_n in 1/19 gen ed = . gen lb = . gen ub = . local j = 1 forvalues i = 0.05(0.05)1 { di `i' estimates restore log4 local lodds = ln(`i'/(1-`i')) nlcom `lodds'/[b2]_b[_cons] + [b3]_b[_cons], post replace ed = _b[_nl_1] in `j' replace lb = _b[_nl_1] - invnormal(.975)*_se[_nl_1] in `j' replace ub = _b[_nl_1] + invnormal(.975)*_se[_nl_1] in `j++' } sum x local min = r(min) local max = r(max) levelsof x twoway rarea lb ub a || /* */ line ed a, /* */ clpattern(solid) /* */ ytitle(ED_100*a) /* */ legend(order( 1 "ci" 2 "ED"))/* */ ymticks(`r(levels)', /* */ tpos(inside) /* */ tlength(*2)) /* */ yline(`min' `max', lpattern(dot)) *----------------------- end example --------------------- (For more on how to use examples I sent to the Statalist, see http://home.fsw.vu.nl/m.buis/stata/exampleFAQ.html ) Hope this helps, Maarten ----------------------------------------- Maarten L. Buis Department of Social Research Methodology Vrije Universiteit Amsterdam Boelelaan 1081 1081 HV Amsterdam The Netherlands visiting address: Buitenveldertselaan 3 (Metropolitan), room Z434 +31 20 5986715 http://home.fsw.vu.nl/m.buis/ ----------------------------------------- ___________________________________________________________ Want ideas for reducing your carbon footprint? Visit Yahoo! For Good http://uk.promotions.yahoo.com/forgood/environment.html * * For searches and help try: * http://www.stata.com/support/faqs/res/findit.html * http://www.stata.com/support/statalist/faq * http://www.ats.ucla.edu/stat/stata/

* * For searches and help try: * http://www.stata.com/support/faqs/res/findit.html * http://www.stata.com/support/statalist/faq * http://www.ats.ucla.edu/stat/stata/

**References**:

- Prev by Date:
**st: estout and likelihood ratio test** - Next by Date:
**Re: st: conditional logistic** - Previous by thread:
**Re: st: reverse prediction - confidence interval for x at given y in nonlinear model** - Next by thread:
**Re: st: reverse prediction - confidence interval for x at given y in nonlinear model** - Index(es):

© Copyright 1996–2024 StataCorp LLC | Terms of use | Privacy | Contact us | What's new | Site index |