How to guess starting value for non-linear regression

Question

I used proc nlmixed in SAS to calculate the beta estimates in non-linear regression model, while I'm not sure how to guess the starting values of the parameters. By using the example from SAS user's Guide, the SAS dataset follows:

data rats;
   input trt $ m x @@;
   if (trt='c') then do;
      x1 = 1;
      x2 = 0;
   end;
   else do;
      x1 = 0;
      x2 = 1;
   end;
   litter = _n_;
   datalines;
c 13 13   c 12 12   c  9  9   c  9  9   c  8  8   c  8  8   c 13 12   c 12 11
c 10  9   c 10  9   c  9  8   c 13 11   c  5  4   c  7  5   c 10  7   c 10  7
t 12 12   t 11 11   t 10 10   t  9  9   t 11 10   t 10  9   t 10  9   t  9  8
t  9  8   t  5  4   t  9  7   t  7  4   t 10  5   t  6  3   t 10  3   t  7  0
;

Assuming I don't know the starting values and thereby set the values as one as following:

proc nlmixed data=rats;
   parms t1=1 t2=1 s1=1 s2=1;
   eta = x1*t1 + x2*t2 + alpha;
   p   = probnorm(eta);
   model x ~ binomial(m,p);
   random alpha ~ normal(0,x1*s1*s1+x2*s2*s2) subject=litter;
   estimate 'gamma2' t2/sqrt(1+s2*s2);
   predict p out=p;
run;

The results from SAS is shown as following:

And then how I can further adjust the starting value based on these results? I have read the method of the grid search and the method of reduced models. Can I just based on the above parameter estimates to adjust the starting value?

Furthermore, I also want to ask that, in some cases, I change the starting value to different number but SAS cannot give me results. It only shows specifications table (like the below pic) but no further results. Is this because my starting value doesn't make sense?

Answer 1

I prefer to avoid the bounds statement whenever possible. At the boundary constraints, the parameters are no longer continuous and derivatives become a problem.

So, what can you do to ensure that the parameters do not violate range restrictions while ensuring that all parameters have derivatives that can be computed at all points of the domain? Reparameterize the model. For standard deviations, employ as parameters log_s1 and log_s2. You can then compute s1 and s2 as exp(log_s1) and exp(log_s2).

Thus, my preferred code (using the same starting values as you specified) would be:

proc nlmixed data=rats gconv=0 absgconv=0 technique=nrridg;
   parms t1=1 t2=1 log_s1=%sysfunc(log(0.05)) log_s2=0;
   s1  = exp(log_s1);
   s2  = exp(log_s2);
   eta = x1*t1 + x2*t2 + alpha;
   p   = probnorm(eta);
   model x ~ binomial(m,p);
   random alpha ~ normal(0,x1*s1+x2*s2) subject=litter;
   estimate 'gamma2' t2/sqrt(1+s2*s2);
   predict p out=p;
run;

The parameters log_s1 and log_s2 have domain from -infinity to infinity, but the values of s1 and s2 have range 0 to infinity. The derivatives of log_s1 and log_s2 can be computed at all values of the parameters.

You will notice, too, that in order to impose the same initial value for s1(=0.05) as you did, I used log_s1=%sysfunc(log(0.05)). Log(0.05) was passed as an argument to the sysfunc function. The sysfunc function is part of SAS scripting capabilities. SAS scripting allows code to be constructed before the code is sent to the SAS compiler. This makes it easy to initialize the value of s1 as 0.05 even though s1 is obtained as exp(log_s1).

Answer 2

In this case you likely don't need to do anything else as all of the gradients are close to zero. Additional iterations would likely not produce different results as the gconv gradient threshold would kick in almost immediately.

You should add in bounds s1 >= 0, s2 >= 0 as those represent non-negative standard deviations and while that doesn't affect the fit of the model in this case because s1*s1 and s2*s2 are then always positive, it is still a good idea. I'd also use the option technique=nrridg. I've just found that option tends to produce smaller gradients.

But an approach to getting starting values for the fixed effects is to fit a model without random effects with a procedure that doesn't require you to set initial values. For your analyses that might be the following:

proc glimmix data=rats;
  class trt;
  model x/m = trt / dist=binomial noint link=probit solution;
run;

This results in

                          Parameter Estimates

                             Standard
Effect    trt    Estimate       Error       DF    t Value    Pr > |t|
trt       c        1.2744      0.1355       30       9.40      <.0001
trt       t        0.7468      0.1153       30       6.47      <.0001

For the variances of the random effects, you might still need to guess.

However...PROC GLIMMIX can also fit the model you've described with PROC NLMIXED.

proc glimmix data=rats method=laplace;
  class trt litter;
  model x/m = trt / dist=binomial noint link=probit solution;
  random intercept / subject=litter group=trt;
run;

Although you'd need to modify your code slightly to estimate the variances rather than the standard deviations to get pretty much the same results:

proc nlmixed data=rats gconv=0 absgconv=0 technique=nrridg;
   parms t1=1 t2=1 s1=0.05 s2=1;
   bounds s1 > 0, s2 > 0;
   eta = x1*t1 + x2*t2 + alpha;
   p   = probnorm(eta);
   model x ~ binomial(m,p);
   random alpha ~ normal(0,x1*s1+x2*s2) subject=litter;
   estimate 'gamma2' t2/sqrt(1+s2*s2);
   predict p out=p;
run;