Problem with nlme() convergence in a Monte Carlo simulation

#1
Hello.

As part of a Monte Carlo simulation, I'm programming some model comparisons on generated data. The underlying relationship in the data is a two-parameter exponential curve (y=a*e^(bx)), and among the models being compared are hyperbolic (y=a/(1+bx)), and quadratic (y=b*(x-a)^2) relationships. The generated datasets follow a within-subjects design, and each subject may have unique values for the two parameters. As such, the models are being fit using nlme().

The problem I'm having is that the hyperbolic and quadratic models fail to converge a majority of the time. The error I usually get for the hyperbolic model is:

Code:
Error in nlme.formula(y ~ a/(1 + b * x), data = theData, fixed = a + b ~  : 
  Step halving factor reduced below minimum in PNLS step
The quadratic model also produces the same type of error, but it will sometimes produce this one instead:
Code:
Error in FUN(newX[, i], ...) : 
  NA/NaN/Inf in foreign function call (arg 1)
Choosing better starting parameter values helps a little, but the problem is still prevalent. I've tried raising the nlmeControl parameters "tol," "pnlsTol," and "minScale," to no avail, and my google-fu hasn't uncovered any other possible solutions.

Here is my code for the model fits:
Code:
hypMod1<-try(nlme(y~a/(1+b*x), data=theData, fixed=a+b~1, random=a+b~1|subj, groups=~subj, start=list(fixed=c(a=ahyp,b=bhyp)), correlation=NULL, control=nlmeControl(maxIter=100, tol=1, pnlsTol=.1, minScale=.1), weights=varPower(), method="REML"))

quadMod1<-try(nlme(y~b*(x-a)^2, data=theData, fixed=a+b~1, random=a+b~1|subj, groups=~subj, start=list(fixed=c(a=aquad,b=bquad)), correlation=NULL, control=nlmeControl(maxIter=100, tol=1, pnlsTol=.1, minScale=.1), weights=varPower(), method="REML"))
(ahyp, bhyp, aquad, and bquad are products of separate estimation functions, similar to a SelfStart approach)

I've attached three example datasets (comma-separated .txt, since I can't do .csv), if that helps. I can't seem to get one that produces the Na/NaN/Inf error for the quadratic model right now, since that one is less common.

Some other info that might be useful: The simulation tests several different ways of sampling the independent variable in the generated datasets. The quadratic model fails more often when a small number (2-4) of discrete IV values are sampled. The hyperbolic model fails more often when there are fewer observations per subject and/or when there is greater between-subject variance in the b parameter in the underlying exponential relationship. Unfortunately, these are all aspects of the simulation that need to remain in place.

I'm pretty desperate at this point. I would be grateful for any help you can offer. Thank you.
 
Last edited:

JesperHP

TS Contributor
#2
When the error is
Code:
Step halving factor reduced below minimum in PNLS step
and you say
I've tried raising the nlmeControl parameters "tol," "pnlsTol," and "minScale,"
I'm thinking you should lower pnlsTol which as I understand it is the minimum in PNLS step?

Still in general you will probably need to sort out some percentage of non-converged modelruns...
 
#3
Thanks for your reply.

Raising pnlsTol was a possible solution I found when googling that error. I don't know a lot about it, but I think it determines how easy it is for the fit to be "good enough" to converge. Nevertheless, I just tried a small run with pnlsTol=.0001 (the default is .001), and it didn't seem to help.

Still in general you will probably need to sort out some percentage of non-converged modelruns...
Oh, definitely. I've actually done simulations like this before (though never with this much trouble), and I can accept non-convergence rates of around 5-7% (perhaps a bit higher if the rate isn't related to my IVs in the simulation). Unfortunately, right now I'm getting non-convergence rates of about 50% for the hyperbolic model and 60% for the quadratic.
 

JesperHP

TS Contributor
#4
I do not know the internals of nlme() but data scaling might be a thing to consider - although perhaps the nlme() call does this automatically - but if it does not then I would consider playing around with the units of the x variable for example using e^(bx) when x is taking on values like 40-60 is not the best idea values around 1=bx is probably preferable....
 
#6
Scaling down the x variable seems to have worked!

I remember reading the reply on that page, but I must have glossed over the OP (and yes, I feel rather foolish for that). Before this, I thought that the scale of the x variable was arbitrary, so this is great to know.

The ideal scale seems to be about one fifth of what I was using (0-20 instead of 0-100). That change fixed everything except the quadratic model's problems when there were 2 or 3 IV levels, but I fixed that by improving my starting parameter value estimating functions. Non-convergence rates for both the hyperbolic and quadratic models are now around 5%.

Thanks so much for your help!

P.S. I saw another thread with a [SOLVED] tag on the title line. Do you know how I can add one of those to this thread?