--- 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/
