[R-sig-ME] Formula and Start-values for simplest "non-linear" mixed model - nlmer (nlme??, lmer ?)

[Ben Bolker]

   Here's a start (I hope) -- new examples for nlmer ...  for the naming
issue, I would consider brute force:  R has a feature that often annoys
me that the names of vectors get carried along ...


> x <- c(a=1,b=2)
> y <- c(a=x[1])
> y
a.a
  1

instead try y<-c(a=unname(x[1]))


> y
a
1

## nonlinear mixed models --- 3-part formulas ---
## 1. basic nonlinear fit. Use stats::SSlogis for its
## implementation of the 3-parameter logistic curve.
## "SS" stands for "self-starting logistic", but the
## "self-starting" part is not currently used by nlmer ... 'start' is
## necessary
startvec <- c(Asym = 200, xmid = 725, scal = 350)
(nm1 <- nlmer(circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym|Tree,
             Orange, start = startvec))
## 2. re-run with "quick and dirty" PIRLS step
(nm1a <- update(nm1, nAGQ = 0L))
## 3. Fit the same model with a user-built function:
## a. Define formula
nform <- ~Asym/(1+exp((xmid-input)/scal))
## b. Use deriv() to construct function:
nfun <- deriv(nform,namevec=c("Asym","xmid","scal"),
              function.arg=c("input","Asym","xmid","scal"))
nm1b <- update(nm1,circumference ~ nfun(age, Asym, xmid, scal)  ~ Asym |
Tree)
## 4. User-built function without using derivs():
##    derivatives could be computed more efficiently
##    by pre-computing components, but these are essentially
##    the gradients as one would derive them by hand
nfun2 <- function(input, Asym, xmid, scal) {
    value <- Asym/(1+exp((xmid-input)/scal))
    grad <- cbind(Asym=1/(1+exp((xmid-input)/scal)),
              xmid=-Asym/(1+exp((xmid-input)/scal))^2*1/scal*
                    exp((xmid-input)/scal),
              scal=-Asym/(1+exp((xmid-input)/scal))^2*
                     -(xmid-input)/scal^2*exp((xmid-input)/scal))
    attr(value,"gradient") <- grad
    value
}
stopifnot(all.equal(attr(nfun(2,1,3,4),"gradient"),
                    attr(nfun(2,1,3,4),"gradient")))
nm1c <- update(nm1,circumference ~ nfun2(age, Asym, xmid, scal)  ~ Asym
| Tree)

On 14-01-16 01:01 PM, Robert Chatfield wrote:
> Hi, Ben.  There's been no response yet to the question about nlmer.   
> Could you start me off with  an answer to the question of how to send
> start values   start.vals=    as a named vector?  
> 
> It's fine to take the answer on the list.
> 
> 
> (NOTE: It does not seem to work to look for the usual  y intercepts by running the problem "backwards"
> as a multiple regression of a single tracer-variable versus many species.type concentrations
> (y[ ij ] == a[ i ] *( + {x^o}[ j ]  
> x[j]  = alpha(i) * y[ij] + x^o[j]  
> since there is an "exact" but trivial solution x^0[j] = x[j]  ... 
> 
> Two questions ... 
> ======================================================================================
> Mostly: (1)
> I could not convince nlmer that the the starting values were vectors ... perhaps this is
> elementary R
> start.vals.1 = c(a=aij, x.nought=x.nought.lin)    (SEE VALUES BELOW)
> #   start.vals.1 = cbind(a=aij, x.nought=x.nought.lin)  # tried this -- does not work
> 
> where this calculation for aij
> len.lon.dat = dim(good.est.merge.long.dat)[1]   #  repeats once for each species.type
> len.lon.dat = dim(good.est.merge.long.dat)[1]   #  repeats once for each species.type 
> 
> aij = vector(mode="numeric",length= len.lon.dat)
> #  The name aij is given since
> for ( i in 1:len.lon.dat ) {
> spec.type.nm = as.character(species.type[i])
> aij[i] = (ranef(try.spec.2.lmer,drop=T)[["species.type"]])[[spec.type.nm]]
> }
> 
> 
>  I've been using c( , ) but I'm wondering
> if nlmer expects the results that occur for a vector of start.vals
> 
> I am wondering if the names are not being passed along correctly ...
>> names(start.vals.1)[1]
> [1] "a1"
>> start.vals.1[1]
>       a1 
> 0.002266363 
>> start.vals.1[1802]
>    a1802 
> -0.02694557 
>> start.vals.1[1803]
> x.nought1 
> -0.3141149 
>> start.vals.1[3604]
> x.nought1802 
> -0.1204595 
> Could this naming be confusing to nlmer?  The start values must be named.
> 
> == from a follow-up post ====  "not all parameter names"  meaning?
> 
>  first.nlmer = nlmer( y ~ a*(x+ x.nought) ~ (a|species.type) + (x.nought|id),
> 
> +      start=start.vals.1,  data=good.est.merge.long.dat  )
> Error in nlmer(y ~ a * (x + x.nought) ~ (a | species.type) + (x.nought |  : 
>   not all parameter names are used in the nonlinear model expression
> 
> where named list is given by ...
>>
>  start.vals.1 = c(a=aij, x.nought=C.nought.lin)
> 
>>
>  start.vals.1[1]
> 
>          a1 
> 0.002266363 
>>
>  start.vals.1[3604]
> 
> x.nought1802 
>   -0.1204595 
> 
> Using single values with a single name also does not work, nor does a construct
> with cbind:
> 
>>
>  first.nlmer = nlmer( y ~ a*(x+ x.nought) ~ (a|species.type) + (x.nought|id),
> 
> +      start=c(a=0.00364901,x.nought=0),  data=good.est.merge.long.dat  )
> Error in nlmer(y ~ a * (x + x.nought) ~ (a | species.type) + (x.nought |  : 
>   gradient attribute of evaluated model must be a numeric matrix
> 
> Bob Chatfield
> 
>