How to specify upper and lower parameter bounds and complex random effects in nlmer models of lme4 when using the bobyqa optimizer

Question

I have a dataset to which I want to fit a nonlinear model with random effects. The dataset involves different lines being observed along time. The total number of lines were split up into batches that were executed on different times in the year. I ran into some issues along the way of analyzing the dataset:

1. how to specify boundaries of parameters using nlmer and the bobyqa optimizer.
2. how to specify complex random effects
3. how to compare nlmer objects with a model without random effects

A simple version of my dataset is as follows:

batch<-c(rep("A",29),rep("B",10),rep("C",10))
line<-c(rep(1:3,9), 1,3,rep(4:5,5),rep(6:7,5))
day<-c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
       2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 
       4L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 7L, 7L, 8L, 8L, 8L, 
       9L, 9L)

result<-c(-2.5336544395728, -2.69179853934892, -2.85649494251061, -4.08634506491338, 
                 -3.57079698958629, -2.62038994824068, -2.69029745619165, -2.18131299587959, 
                 -2.1028751114459, -2.56553968316024, -2.55450633017557, -2.43072209061048, 
                 -2.42496349148255, -2.52850292795008, -1.09958807849945, -1.49448455383069, 
                 -0.461525929110392, -0.298569396331159, -0.520372425046126, -0.676393841451061, 
                 -0.930448741799686, -0.414050789171074, -0.0915696466880981, 
                 -0.239509444743891, -0.319036274966057, -0.189385981406834, -0.376015368786241, 
                 -0.269728570922294, -0.260869642513491, -0.206260420960064, -0.790169432375232, 
                 -0.0573210164472325, -0.202013642441365, -0.0853200223702248, 
                 -0.13422881481997, 0.0831839881028635, -0.0288333371533044, 0.124233139837959, 
                 -0.16906818823674, -0.299957519185662, -0.085547531863026, 0.00929447912629702, 
                 -0.117359669415726, -0.0764263122102468, -0.00718772329252618, 
                 0.0110076995240481, -0.0304444368953004, 0.0586926009563272, 
                 -0.0457574905606751)

data <- data.frame(day, line, batch, result)
data$line<-as.factor(data$line)
data$batch<-as.factor(data$batch)

The nonlinear model without random effects fits just fine and is defined as:

nlsfit=nls(result~ z-(p0*(z-Za))/(p0+(1-p0)*(1/(1+s))^day), 
       data=data, start=list(z=-8,p0=0.0001,Za=0,s=10), algorithm="port")

And gives following output:

Nonlinear regression model model: result ~ z - (p0 * (z - Za))/(p0 + (1 - p0) * (1/(1 + s))^day) data: data z p0 Za s -3.10281 0.02946 -0.12286 9.49270 residual sum-of-squares: 4.158

Algorithm "port", convergence message: relative convergence (4)

When incorporating mixed effects caused by batch and line factors, I realized that the nlme() function from the nlme package has issues with a more complex structure of the mixed models (or I failed to implement it correctly). However, a fairly simple implementation of a random effect works ànd I can specify the upper and lower limits easily:

upper_bounds<-c(2, 1, 2, 50)
lower_bounds<-c(-8, 0.00000001,-8,0)
f1 <- result ~ z-(p0*(z-Za))/(p0+(1-p0)*(1/(1+s))^day)

nlmefit=nlme(f1,
             data = data,
             fixed = z + p0 + Za + s ~ 1, 
             random = z ~ 1|batch,
             start = coef(nlsfit),
             control = nlmeControl(opt = "nlminb",                               
                                   lower = lower_bounds, 
                                   upper = upper_bounds, 
                                   maxIter = 1000, 
                                   msMaxIter = 1000))

However, the logically expected structure of random effects is much more complex. Specifically, more parameters than simply "z" are likely affected by batch and/or line effects but already adding one additional parameter to the random effect, results in convergence/singularity errors. Luckily, the nlmer() function of lme4 does allow for more complex random effects to be specified. Furthermore, when using bobyqa as optimizer, my nlmer model has no convergence issues. E.g.:

#defining the function needed for nlmer()
nform <- ~ z-(p0*(z-Za))/(p0+(1-p0)*(1/(1+s))^day)
nfun <- deriv(nform, namevec = c("z","p0","Za","s"),
              function.arg = c("day", "z","p0","Za","s"))

nlmerfit = nlmer(log10perfract ~ nfun(day, z, p0, Za, s) ~ 
                   (z+s+Za|batch),
                 data = data,
                 start= coef(nlsfit),
                 control= nlmerControl(optimizer = "bobyqa")

1. However, specifying upper and lower limits does not work as before with nlme. resulting in:

Error in nlmerControl(optimizer = "bobyqa", lower = lower_bounds, upper = upper_bounds) : unused arguments (lower = lower_bounds, upper = upper_bounds)

When specifying these bounds in an optCtrl argument as a list, R returns that my starting values violate the bounds (which they do not?):

nlmerfit = nlmer(log10perfract ~ nfun(day, z, p0, Za, s) ~ 
                   (z+s+Za|batch),
                 data = data,
                 start= coef(nlsfit),
                 control= nlmerControl(optimizer = "bobyqa", 
                                       optCtrl = list(lower = lower_bounds,
                                                      upper = upper_bounds)
                                       )
                  )

Error in (function (par, fn, lower = -Inf, upper = Inf, control = list(), : Starting values violate bounds

I need these bounds to be working as with nlme as my real data is a bit more complex (containing different groups of data for which the bounds are needed to allow a fit).

2. How to specify more complex random effects like nested effects? For example, both s, Za and p0 are expected to show some random variation between individual lines. Similarly, z is expected to vary between lines, but much more so between lines of different batches (nested?). The batch effect on s, Za and p0 are less logic from a biological/evolutionary point of view but, looking at my data, seems to be present nevertheless. What is the right syntax to implement these effects?

3. How can I compare nlmer models with nlme or nls models? using anova(object1,object2) I can compare different nlmer objects or compare the nlme model with simple random effects with the nls fit without such effects. When I try to compare the nlmer model with either nlme or nls models, I get the following error:

anova(nlmerfit,nls)

Error in .rowNamesDF<-(x, value = value) : duplicate 'row.names' are not allowed In addition: Warning messages: 1: In anova.merMod(nlmerfit, nls) : additional arguments ignored 2: non-unique values when setting 'row.names':

Any help is much welcome!