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

["Rosy Reynolds" <rr@dandr.demon.co.uk>]

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" <maartenbuis@yahoo.co.uk>
To: <statalist@hsphsun2.harvard.edu>
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 <rr@dandr.demon.co.uk> 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/