[R-sig-ME] Formula and Start-values Any help at all?… Self-contained example

[Ben Bolker]

Robert Chatfield <chatfield at ...> writes:

>  Dear list, I still have no idea how to specify or give start values
> for vectors in a formula for nlmer() The technique works fine for an
> approximate method I use in lmer() (That technique depends heavily
> on some scaling of variables.)  Here is a self-contained example
> using some mock data whose origin is described.

  I don't know if I got sensible answers or not, but I seem to have
succeeded in fitting a model.  The key changes were: (1) starting
parameters are *single* values (these are starting values for the
mean/central values of the parameters, not of the random effects);
(2) leaving "xlong" out of namevec in the deriv() call, because it's
a predictor variable, not a parameter (you don't want its gradient)
                                             
Slight modifications:

## Example of emissions estimation

set.seed(101)
## Define a "true" background concentration ... 100 instnaces
x0 = rnorm(100, 20, 5)
## Define a perturbation of the key tracer -- log normal
xd = rlnorm( 100, log(8), 0.3)
## Define the measurable, x
x = x0 + xd
## Define some "emission factors" linearly relating several dependent variables 
## as they are related to the key tracer (above background)
## 10 different dependent variables or "species"
a = 0.1*rexp(10,rate=0.6)
## Define the pure response variable (yun ... "y unscaled")
yun = t( outer( a,xd ) )
## Add some noise in both the response, both w.r.t. the dependent species 
##    and also the measurement instance
## for (i in 1:10 ) yun[,i] = yun[,i]*(1+rnorm(100,0.05,0.1))
## for (j in 1:100 ) yun[j,] = yun[j,]*(1+rnorm(10,0.05,0.1))
yun <- sweep(yun,1,rnorm(nrow(yun),0.05,0.1),"*")
yun <- sweep(yun,2,rnorm(ncol(yun),0.05,0.1),"*")
##
y = yun    
##  dimension a scaled response and fill it with scaled values
##    (still all positive)
ysc = vector(mode="numeric", length=10)
asc = vector(mode="numeric", length=10)
for ( i in 1:10 ) { 
    ytemp =  scale(yun[,i], center=FALSE) # get ytemp to define scaling factor
    y[,i] =  ytemp   # scale
    ysc[i] = attributes(ytemp)$`scaled:scale`
    asc[i] = 1/ysc[i]    # useful version of a to relate scaled quantities
}
##
## Make a "longi" list of the tracer-species, repeating once for every response
## Note:  normally I use reshape() to do this.  See: ?reshape
xlong = rep(x,times=10) 
ylong = as.vector(y) # collapse the matrix into a similar "long" vector
## Make a factor variable describing each measurement 
## instance and each dependent variable 
id = as.factor(rep(paste("sample",as.character(1:100),sep=""),times=10))
species = as.factor(rep(paste("species",as.character(1:10),sep=""),times=100))
##
## Create the data-frame for nlmer
trial.dat = data.frame(xlong=xlong,ylong=ylong,id=id,species=species)
##
xinterform <- ~( a * xlong+ x0.est)
xinter <-  deriv(xinterform, namevec=c("a","x0.est"),
                 function.arg=c("a","xlong","x0.est"))
##
library(lme4)
start.vec <- c(a=0.16,x0.est=20)
##
test.nlmer = nlmer( ylong ~ 
    xinter(a,xlong,x0.est)  ~ (a|species) + (x0.est |id),
    start=start.vec,  data=trial.dat  )