[R-sig-ME] Formula and Start-values Any help at all?… Self-contained example
Ben Bolker
bbolker at gmail.com
Sat Jan 25 20:45:05 CET 2014
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 )
More information about the R-sig-mixed-models
mailing list