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

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

From   "Rosy Reynolds" <>
To   "statalist" <>
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
BSAC Resistance Surveillance Co-ordinator

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

----- Original Message ----- From: "Maarten buis" <>
To: <>
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 <> 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

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:

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

Hope this helps,

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

Want ideas for reducing your carbon footprint? Visit Yahoo! For Good
*   For searches and help try:

*   For searches and help try:

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