David Airey was asking about the ability of Stata to implement nonlinear curve fitting for
selected nonlinear functions of interest to pharmacologists. I had posted that Stata has
-nl- and -ml- that would be helpful in this regard, and that -ml- has certain advantages,
for example, the -cluster()- option. He replied,
-------------------------begin excerpt from post---------------------------------
Thanks for this information. I'll likely give one of these equations a
whirl to compare results. Others in the lab uses Motulsky's GraphPad
prism because it's easy to use and the documentation is well done if
full of typos. Be nice to have extended capabilities present in Stata
that are not present in Prism.
From what I can tell, Prism 3.0 has several built in (and modifiable)
classic equations, such as:
1. one site binding (hyperbola)
Y = Bmax * X
--------
Kd + X
two site binding
Y = Bmax1 * X Bmax2 * X
--------- + ---------
Kd1 + X Kd2 + X
2. sigmoidal dose-response
(top - bottom)
Y = bottom + --------------------
1 + 10^(logEC50 - X)
3. sigmoidal dose-response (variable slope)
(top - bottom)
Y = bottom + ------------------------------
1 + 10^(logEC50 - X)*hillslope
4. one site competition
(top - bottom)
Y = bottom + --------------------
1 + 10^(X - logEC50)
5. two site competition
fraction_1 1 - fraction_1
Y = bottom + (top - bottom) *( --------------------- + -------------------- )
1 + 10^(X-logEC50_1) 1 + 10^(X-logEC50_2)
and a few others:
6. boltzmann sigmoid
7. one phase exponential decay
8. two phase exponential decay
9. one phase exponential association
10. two phase exponential association
11. exponential growth
12. power series
13. polynomial equations
-------------------------end excerpt from post----------------------------------
A couple of items from the list, for example, the exponential decay and association
functions (Nos. 7 through 11), might already be programmed for -nl-; there is a suite of
exponential functions for -nl- that ship with Stata. It looks as if all of the functions
through Number 11 would be easy to implement with -nl- if they're not already
available. The last two, I believe, are currently implemented or served in various ways
in Stata, for example, in -boxcox-, -fracpoly-, -boxtid- and -bspline-.
In case it helps to illustrate the use of -ml- for these kinds of problems, I've
programmed the first two from the list (one-site and two-site Langmuir adsorption
isotherms) and attached them below in a do-file. The illustrations below use artificial
datasets for a problem set exercise in a biochemistry course taught by Richard Neubig
at the University of Michigan. They can be found as the first two hits from Google
searching using keywords, *neubig* and *550*. (The URLs are each too long to include
either on a single line here.)
If there is an interest in implementing nonlinear regression items from the list above
using -ml-, good references for the exercise are:
1. help for -ml-,
2. Section 3.4 "Nonlinear specifications" in W. Gould & W. Sribney, _Maximum
Likelihood Estimation in Stata_ (College Station, Texas: Stata Press, 1999), pp. 44-45
[note that there's a typographical error in the code shown on Page 45--a right
parenthesis is missing in each case], and
3. Chapter 13. "Maximum likelihood estimation" in S. Rabe-Hesketh & B. Everitt, _A
Handbook of Statistical Analyses usng Stata_, Second Edition. (Boca Raton, Florida:
Chapman & Hall/CRC, 2000), esp. pp. 181-82.
Joseph Coveney
--------------------------------------------------------------------------------
/* estimating receptor binding parameters at equilibrium;
one binding site
variables are fre (free ligand), tbn (total
binding), nsb (nonspecific binding) and nbi
(net or specific binding) */
clear
set more off
input long fre int tbn int nsb
5000 325 25
10000 550 50
20000 900 100
50000 1300 200
100000 1700 400
200000 2200 800
500000 3200 1600
end
/* conversion from cpm to femtomoles */
foreach var of varlist _all {
replace `var' = `var' / 0.331 / 2.22 / 136
}
generate float nbi = tbn - nsb
*
*
program define mlonesite
args lnf theta1 theta2 theta3
tempvar res
quietly generate double `res' = $ML_y1 - /*
*/ fre * `theta1'/ (`theta2' + fre)
quietly replace `lnf' = -0.5 * ln(2*_pi) - /*
*/ ln(`theta3') - 0.5 * `res'^2 / `theta3'^2
end
*
*
ml model lf mlonesite (Bmax:nbi=) (Kd:) (sigma:)
ml check
ml search
ml maximize, difficult
*
*
*
/* two binding sites
variables are fre (free ligand) and nbi
(net or specific binding) */
clear
input float fre float nbi
0.25 12.4
0.50 21.4
0.75 28.4
1.0 34.1
1.5 43
2.0 50
3.0 61
5.0 75
10 96
15 107
20 114
25 120
end
aformat _all
clist
*
*
program define mltwosite
args lnf theta1 theta2 theta3 theta4 theta5
tempvar res
quietly generate double `res' = /*
*/ $ML_y1 - ( fre * `theta1'/ (`theta2' + fre) + /*
*/ fre * `theta3' / (`theta4' + fre) )
quietly replace `lnf' = -0.5 * ln(2*_pi) - /*
*/ ln(`theta5') - 0.5 * `res'^2 / `theta5'^2
end
*
*
ml model lf mltwosite (Bmax1:nbi=) (Kd1:) (Bmax2:) (Kd2:) (sigma:)
ml check
ml search
ml maximize, difficult
exit
--------------------------------------------------------------------------------
*
* 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/
