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

[Robert Chatfield]

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 always get one of these two errors when trying to specify the nlmer() model
Error in nlmer(ylong ~ xinter(a, xlong, x0.est) ~ (a | species) + (x0.est |  : 
  not all parameter names are used in the nonlinear model expression
or
Error: s > 0 is not TRUE

Can  I get any help on formula specification or start values?
I have coded in a function with derivatives using deriv() .


> # Example of emissions estimation
> #
> # 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))
> #
> 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=F) # 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)
> #
> # Define a formula ... trying to force a vector definition 
> ##  xinterform <- ~a*(xlong+ x0.est)  # previous try
> #
> xinterform <- ~( a * xlong+ x0.est)
> # Explicitly declaring these as vectors gives an error
> # xinterform <- ~( a=vector(mode="numeric",length=1000) * xlong=vector(mode="numeric",length=1000) + x0.est=vector(mode="numeric",length=1000) )
> # Doesn't work:
> # Error in deriv.formula(xinterform, namevec = c("a", "xlong", "x0.est"),  : 
> #  invalid expression in 'FindSubexprs'
> 
> xinter <-  deriv(xinterform, namevec=c("a","xlong","x0.est"),
+                     function.arg=c("a","xlong","x0.est"))
> #
> start.vec = c(a=rep(0.16,times=1000),x0.est=rep(20,times=1000))
> #
> test.nlmer = nlmer( ylong ~ xinter(a,xlong,x0.est)  ~ (a|species) + (x0.est |id),
+      start=start.vec,  data=trial.dat  )
Error in nlmer(ylong ~ xinter(a, xlong, x0.est) ~ (a | species) + (x0.est |  : 
  not all parameter names are used in the nonlinear model expression
> #
> start.vec = unname(start.vec)
> test.nlmer = nlmer( ylong ~ xinter(a,xlong,x0.est)  ~ (a|species) + (x0.est |id),
+      start=start.vec,  data=trial.dat  )
Error: s > 0 is not TRUE

Notes -- how R defines start.vec...


> start.vec = unname(start.vec)
> start.vec[c(1,1000, 1001,2000)]
[1]  0.16  0.16 20.00 20.00

> start.vec = c(a=rep(0.16,times=1000),x0.est=rep(20,times=1000))
> start.vec[c(1,1000, 1001,2000)]
        a1      a1000    x0.est1 x0.est1000 
      0.16       0.16      20.00      20.00