[R] nlsList {nlme} - control arguments problem

[Rick DeShon]

Hi All.

I'd like to send some control arguments to the nls function when
performing a nlsList analysis.

I'm fitting a power model to some grouped data and would like to
impose lower bounds on the estimates using the "port" algorithm.
Obtaining the lower bound constraint works fine with a direct call to
nls for a single level of the grouping variable.  However, the bounds
aren't imposed when performing the nlsList analysis across the levels
of the grouping variable.  Any idea why this isn't working?

# Generate example data ##########

trial  <- 1:100
result <- list()

for (i in 1:3) {
 min <- rnorm(max(trial),250,5)
 dif <- rnorm(max(trial),500,5)
 p   <- rnorm(max(trial),-.5,.1)
 e   <- rnorm(max(trial),mean=0,sd=5)
 y   <- min + dif*(trial)^p + e
 result[[i]] <- data.frame(y,trial,id=i)
}
newdf   <-do.call('rbind',result)
df.gr     <- groupedData( y ~ trial | id, data=newdf)


### Single unit analysis
........................................................................
### The boundary condition on the dplt parameter is enforced! ..........

df.one <- subset(df.gr,id==1)
nls(y~SSpowrDplt(trial,min,dplt,dif,p),data=df.one,algorithm="port",lower=c(0.0,0.0,0.0,-10))

...... example output.......
>Nonlinear regression model
>  model:  y ~ SSpowrDplt(trial, min, dplt, dif, p)
>  data:  df.one
>    min    dplt     dif       p
>247.052   0.000 491.965  -0.462
> residual sum-of-squares: 234322
> Algorithm "port", convergence message: relative convergence (4)
##################################

### nlsList analysis
........................................................................................
### Boundary condition on dplt is not enforced
.............................................

Lfit.nls    <- nlsList(y~SSpowrDplt(trial,min,dplt,dif,p),data=df.gr,control=list(algorithm='port',lower=c(0.0,0.0,0.0,-10),maxiter=100))

...... example output.......
>Call:
>  Model: y ~ SSpowrDplt(trial, min, dplt, dif, p) | id
>   Data: df.gr

>Coefficients:
>   min
>  Estimate Std. Error   t value     Pr(>|t|)
>1 276.2354   16.16609 17.087337 1.086442e-44
>2 257.0127   20.43564 12.576694 3.390686e-30
>3 206.4017   29.01315  7.114075 7.354863e-12
>   dplt
>     Estimate Std. Error    t value  Pr(>|t|)
>1 -0.06579982 0.03848086 -1.7099365 0.0951222
>2 -0.01694362 0.04161933 -0.4071093 0.6786473
>3  0.08981518 0.04636532  1.9371199 0.0528957
>   dif
>  Estimate Std. Error  t value     Pr(>|t|)
>1 477.5049   21.89002 21.81382 6.679439e-62
>2 488.7171   22.11908 22.09482 4.466288e-66
>3 552.7105   25.04206 22.07129 9.215344e-65
>   p
>    Estimate Std. Error   t value     Pr(>|t|)
>1 -0.5455936 0.06262040 -8.712713 7.615265e-16
>2 -0.4839114 0.06074282 -7.966560 1.307734e-14
>3 -0.4059903 0.05455864 -7.441355 9.297527e-13

>Residual standard error: 27.43384 on 888 degrees of freedom
#####################################################

If you look at the structure of Lfit.nls, it looks like the control
arguments are passed correctly.
str(Lfit.nls)

List of 3
 $ 1:List of 6
 ..$ control    :List of 7
 .. ..$ maxiter  : num 100
 .. ..$ tol      : num 1e-05
 .. ..$ minFactor: num 0.000977
 .. ..$ printEval: logi FALSE
 .. ..$ warnOnly : logi FALSE
 .. ..$ algorithm: chr "port"
 .. ..$ lower    : num [1:4] 0 0 0 -10



If it helps, here's the selfStart function that I'm using....
powrDpltInit <-
function(mCall, LHS, data) {
  xy     <- sortedXyData(mCall[["x"]],LHS,data)
  min.s  <- min(y)
  dif.s  <- max(y)-min(y)
  dplt.s <- 0.5
  p.s    <- -.20
  value  <- c(min.s, dplt.s, dif.s, p.s)
  names(value) <- mCall[c("min","dplt","dif","p")]
  value
}

SSpowrDplt<-selfStart(~min + dplt*x + dif*x^p,initial=powrDpltInit,
parameters=c("min","dplt","dif","p"))



Thanks for your help!

Rick DeShon