[R] nls, reach limit bounds

Wed Jul 15 17:49:51 CEST 2009

... and getting an answer is no assurance that the answer is meaningful. In
cases like this (which arise frequently because of insistence on using the
accepted mechanistic model paradigm even in the absence of informative data
to estimate it), the confidence intervals for (correlated) parameters will
usually be so wide as to be useless. That is, for practical purposes, the
model is nonidentifiable. This can easily be seen by making small random
perturbations in the data and watching how the parameter estimates change --
i.e. performing a sensitivity analysis. Incidentally, the predictions will
typically be fine, so the standard scientific practice of graphing the data
with an overlaid smooth curve as a "check" on the validity of the estimated
parameters is nonsense. 

One should not get so lost among the trees of statistical niceties that one
loses sight of the scientific forest: the model can't be adequately fit by
the data. Either get more data or choose a more appropriate model.

Cheers,

Bert Gunter
Genentech Nonclinical Biostatistics

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Ravi Varadhan
Sent: Tuesday, July 14, 2009 4:35 PM
To: 'UyenThao Nguyen'; 'spencerg'
Cc: r-help at r-project.org
Subject: Re: [R] nls, reach limit bounds

I took a quick look at "drc"package and the "drm" function.  The drm()
function uses optim ("BFGS" method).  So, that is one diffference.  However,
without looking at your code on how you used drm(), I cannot tell further.  

The fact that you got an answer using optim() does not necessarily mean that
everything is swell.  Did you check the Hessian to see if it is
positive-definite?

You might also want to try the nls.lm() function in the "minpack.lm"
package.

Ravi.

----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml

----------------------------------------------------------------------------
--------

-----Original Message-----
From: UyenThao Nguyen [mailto:unguyen at tethysbio.com] 
Sent: Tuesday, July 14, 2009 7:07 PM
To: Ravi Varadhan; 'spencerg'
Cc: r-help at r-project.org
Subject: RE: [R] nls, reach limit bounds

Hi Ravi and Spencer,

Thank you very much for your help.

I did plot the data, and saw that the data didn't seem to have an inflection
point. Yes, the data contained 6 points of duplicates, which the 4 P
logistic regression is appropriate to use. 

I tried the dose response model (drm in drc package), this model works
without any problem. Do you know if the drm has different tolerance in
convergent or something else? Let's say if I choose drm to fit the data, Can
I get the parameters in the same way nls gives me? I tested a converged data
set on both drm and nls, and I can't see any relationship between their
parameters although the fits were similar.

Thank you.
Uyen

-----Original Message-----
From: Ravi Varadhan [mailto:RVaradhan at jhmi.edu]
Sent: Monday, July 13, 2009 3:32 PM
To: 'spencerg'; UyenThao Nguyen
Cc: r-help at r-project.org
Subject: RE: [R] nls, reach limit bounds

Hi Uyen, 

You do not have enough information to estimate 4 parameters in your
nonlinear model.  Even though you have 12 data points, only 6 are
contributing independent information (you essentially have 6 replicate
points).  If you plot the fittted function you will see that it fits your
data really well.  However, you will also see that the fitted curve is far
from reaching the asymptote.  An implication of this is that you cannot
reliably estimate b0 and b1 separately.  So, you need more data, especially
for larger x values in the asymptotic region.

Ravi.

----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml

----------------------------------------------------------------------------
--------

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of spencerg
Sent: Saturday, July 11, 2009 9:50 PM
To: UyenThao Nguyen
Cc: r-help at r-project.org
Subject: Re: [R] nls, reach limit bounds

      Have you plotted the data?  There are two standard sources of
non-convergence problems like this:  First, there may not be enough
information in your data to estimate all four parameters.  Second, if that's
not the case, then  your starting values are not appropriate. 

      I routinely use "optim" or "nlminb" to find a sensible solution,
because these general purpose optimizers are more likely to converge than
"nls".  To do this, I write a function with a name like "SSElogistic" to
compute the sum of squares of residuals for your model.  I like to use
"optim" with "hessian=TRUE".  Then I compute "eigen(fit$hessian,
symmetric=TRUE)", with "fit" = the object returned by "optim".  If the
smallest eigenvalue is negative, it says that optim found a saddle point.
If the smallest eigenvalue is less than, e.g.,
1e-8 times the largest, it says that the smallest eigenvector is very poorly
estimated.  Round those numbers off grossly to 1 significant digit, and they
will likely suggest which parameter you can fix and drop from the model. 

      Hope this helps. 
      Spencer Graves

UyenThao Nguyen wrote:
> Hi,
>
> I am trying to fit a 4p logistic to this data, using nls function. The
function didn't freely converge; however, it converged if I put a lower and
an upper bound (in algorithm port). Also, the b1.A parameter always takes
value of the upper bound, which is very strange. Has anyone experienced
about non-convergent of nls and how to deal with this kind of problem?
>
> Thank you very much.
>
>
>
> ########################################################################3
>            y           x
> 1  0.8924619 -0.31875876
> 2  1.1814749 -0.21467016
> 3  1.6148266  0.06069784
> 4  2.2091363  0.54032947
> 5  2.7019079  1.04921802
> 6  3.0679585  1.60745502
> 9  0.9436973 -0.31875876
> 10 1.2201133 -0.21467016
> 11 1.6470043  0.06069784
> 12 2.2090048  0.54032947
> 13 2.6864857  1.04921802
> 14 3.0673523  1.60745502
>
> new.cont=nls.control(maxiter = 10000, tol = 1e-05, minFactor = 1e-08,
>             printEval = FALSE, warnOnly = FALSE)
>
>
> b0.A=.9*min(DAT$y)
> b1.A=1.1*max(DAT$y)-b0.A
> b2.A=-1*mean(DAT$x)
> b3.A=1
>
>
> b0.A
> b1.A
> b2.A
> b3.A
>
> nls.mdl.A=nls(y~b0 + b1/(1+exp(-b2-b3*x)),data=DAT,start = 
> list(b0=b0.A, b1=b1.A, b2=b2.A, b3=b3.A), lower=-10, upper=10,
> algorithm="port",trace=T,control=new.cont)
>
> ##################################
>
>
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.