[R] nls() and nls2() behavior?

ivo welch ivowel at gmail.com
Tue May 11 16:46:26 CEST 2010


first, apologies for so many posts yesterday and today.  I am
wrestling with nls() and nls2().  I have tried to whittle it down to a
simple example that still has my problem, yet can be cut-and-pasted
into R.  here it is:

library(nls2)

options(digits=12);

y= c(0.4334,0.3200,0.5848,0.6214,0.3890,0.5233,0.4753,0.2104,0.3240,0.2827,0.3847,0.5571,0.5432,0.1326,0.3481)
x= c(0.3521,0.4334,0.3200,0.5848,0.6214,0.3890,0.5233,0.1379,0.2104,0.3240,0.3404,0.3847,0.5571,0.5432,0.1326)
dummy=c(rep( "m1", 7), rep( "m2", 3), rep( "m3", 5))

d= data.frame( y=y, x=x, dummy=dummy );

print(lm( y ~ as.factor(dummy) + x, data=d ))  ## a simple fixed-effects model


objective.xs <- function( b, m1, m2, m3, x, ydebug ) {
  tgts.mean= c(rep( m1, 7), rep( m2, 3), rep( m3, 5));
  yhat= tgts.mean + b*x
  #debug
    cat("Parameters: b=", b, " m123=", m1, m2, m3, " leading to ",
sum( (ydebug-yhat)^2 ), "\n")
  return( yhat );
}

sameformula= (y ~ objective.xs(b, m1, m2, m3, x, y))

cat("\nNLS2 Function --- what do you do?:\n");
print(nls2( sameformula, data=d, control= nls.control(tol=1e-12),
algorithm="brute-force",
            start=list( b=0.3, m1=0.5, m2=0.5, m3=0.5, trace=TRUE ) ) );

cat("\n---but with this b starting value, this is so much better (?!):\n")
print(nls2( sameformula, data=d, control= nls.control(tol=1e-12),
algorithm="brute-force",
           start=list( b=0.2, m1=0.5, m2=0.5, m3=0.5 ) ));

cat("\nNLS Function:\n");

print(nls( sameformula,  data=d, control= nls.control(tol=1e-12),
          start=list( b=0.3, m1=0.5, m2=0.5, m3=0.5, trace=TRUE ) ));


for some reason, nls2() does not start wandering off into a good
direction; and nls() has no direction (Error in nlsModel(formula, mf,
start, wts) : singular gradient matrix at initial parameter
estimates).  I am probably doing something really silly, but I have
been staring at this example (not just the whittled down, but the
original one from which it was derived) for hours now, and I cannot
figure out how I have miscoded this.  I hope the problem is obvious to
regular users...thanks for any helpful eyes here.

regards,

/iaw
----
Ivo Welch (ivo.welch at brown.edu, ivo.welch at gmail.com)



More information about the R-help mailing list