How to fit a nls model with mixed effects

Question

I want to fit a linear-plateau model with random effects. I found a way to fit the function with nls(), but I don't know how to include random effects. Here is what i have so far:

#create data
x=c(1:6,1:6)
y=c(10,21,27,35,33,35,9,20,28,32,30,31)
z=c("A","A","A","A","A","A","B","B","B","B","B","B")
df<-data.frame(x,y,z)

#create linear-plateau function
lp=function(x, a, b, c){
  ifelse(x > c, a + b * c, a + b * x)
  }

#fit the model without random effects
p10=nls(y ~ lp(x, a, b, c), data = df, start = list(a = 0, b = 15, c = 4))


plot(y~x)
lines(x=c(0, coef(p10)["c"],max(df$x)), 
      y=c(coef(p10)["a"],
          (coef(p10)["a"] + coef(p10)["b"] * coef(p10)["c"]),
          (coef(p10)["a"] + coef(p10)["b"] * coef(p10)["c"])),lty=2)

What I want to do is to include z as a random effect, since all data collected from the same z level are not independent. I know how to model mixed effects with nlmer function from lme4 package, but I don't know how to fit the linear-plateau model with it.

Answer 1

You can do this with the nlme package, but the data you've given us aren't going to be sufficient to succeed in fitting a random-effects model.

Start by fitting a gnls() (generalized non-lin least squares) model, which allows for fixed-effect differences among groups:

library(nlme)
p20 = gnls(y ~ lp(x, a, b, c),
           params= list(a+b~z, c~1),
           data = df,
           start = list(a = c(0,0), b=c(15,15), c=4))

(I initially tried params = list(a+b+c~z), with appropriate changes to start, but the fit didn't succeed. It might be possible to tweak control parameters to make that model work ...)

Now as a random effects model. This doesn't succeed - you would almost definitely need to have more than two groups - but it should give you the idea.

p30 = nlme(y ~ lp(x, a, b, c),
           random = a+b~1|z,
           fixed = a+b+c ~ 1,
           data = df,
           start = c(a=0, b=15, c=4)
           )

Doing this with nlmer is a little fussier as you have to define a function that returns the gradients as well as the value of the objective function.