[R] ask for help on nonlinear fitting
Dieter Menne
dieter.menne at menne-biomed.de
Sat Mar 8 10:13:37 CET 2008
Ruqiang Liang <ruqiang.liang <at> gmail.com> writes:
>
> I have a table like the following. I want to fit Cm to Vm like this:
> Cm ~Cl+Q1*b1*38.67*exp(-b1*(Vm-Vp1)*0.03867)/(1+exp(-b1*(Vm-Vp1)*0.03867))^2+
Q2*b2*38.67*exp(-b2*(Vm-Vp2)*0.03867)/(1+exp(-b2*(Vm-Vp2)*0.03867))^2
>
> I use nls, with start=list(Q1=2e-3, b1=1, Vp1=-25, Q2=3e-3, b2=1,
> Vp2=200). But I always get 'singular gradient' error like this. But
> in SigmaPlot I can get the result. How can I get with R.
I remember a similar case in gastric emptying fitting where my colleague tried
to show me that SigmaPlot is superior to R/nls. When he looked at the standard
deviations of the coefficients, giving estimated gastric empty times of 40
minutes plus minus 800 minutes, I was reminded of Douglas Bates' philosophy:
better nothing than wrong. (Even if I would love to have estimated p-values back
in lmer fits.)
This said, there are two solutions: first, simplify your expression to see what
you are doing. The first part of the expression can be written as
R1*exp(-c1*(Vm-Vp1))/(1+exp(-c1*(Vm-Vp1))^2
where you are free to compute Q1 and b1 later. You will also note that your fit
cannot have a unique solution because of the symmetry in the second term.
As a next step, try to fit
R1*exp(-exp(d1)*(Vm-Vp1))/(1+exp(-exp(d1)*(Vm-Vp1))^2
This gives a working solution in most cases, effectively forcing a positive
c1=exp(d1).
However, even that may fail in degenerate cases. If it is not a single fit, but
from a planned pharmacology related experiment, you should try to fit the whole
experiment with all repeats instead of single curves. This often gives excellent
results even when some fits are disastrous. Check package nlme, the book by
Pinheiro/Bates, and some related examples on http://www.menne-biomed.de/gastempt
Dieter
More information about the R-help
mailing list