# Non-linear mixed effects regression in R

Surprisingly, I was unable to find an answer to the following question using Google:

I have some biological data from several individuals that show a roughly sigmoid growth behaviour in time. Thus, I wish to model it using a standard logistic growth

P(t) = k*p0*exp(r*t) / (k+p0*(exp(r*t)-1))


with p0 being the starting value at t=0, k being the asymptotic limit at t->infinity and r being the growth speed. As far as I can see, I can easily model this using nls (lack of understanding on my part: why can I not model something similar using standard logit regression by scaling time and data? EDIT: Thanks Nick, apparently people do it e.g. for proportions, but rarely http://www.stata-journal.com/article.html?article=st0147 . Next question on this tangent would be if the model can possibly handle outliers >1).

Now I wish to allow some fixed (mainly categorical) and some random (an individual ID and possibly also a study ID) effects on the three parameters k, p0 and r. Is nlme the best way of doing this? The SSlogis model seems sensible for what I am trying to do, is that correct? Is either of the following a sensible model to begin with? I cannot seem to get the starting values right and update() only seems to work for random effects, not fixed ones - any hints?

nlme(y ~ k*p0*exp(r*t) / (k+p0*(exp(r*t)-1)), ## not working at all (bad numerical properties?)
data = data,
fixed = k + p0 + r ~ var1 + var2,
random = k + p0 + r ~ 1|UID,
start = c(p0=1, k=100, r=1))

nlme(y ~ SSlogis(t, Asym, xmid, scal), ## not working, as start= is inappropriate
data = data,
fixed = Asym + xmid + scal ~ var1 + var2, ## works fine with ~ 1
random = Asym + xmid + scal ~ 1|UID,
start = getInitial(y ~ SSlogis(Dauer, Asym, xmid, scal), data = data))


As I am new to non-linear mixed models in particular and non-linear models in general, I would appreciate some reading recommendations or links to tutorials / FAQs with newbie questions.

• If you regard k as known, you could scale your populations by P/k. If k is something to be estimated, that alone means that your problem is not standard logit regression. – Nick Cox Aug 11 '13 at 8:57
• Thank you Nick. Yes, in the end I believe that k needs to be estimated and included in the regression. My interest in using logit regression was purely academic. I thought this might be a nice model to start with before going to non-linear modelling, but I was unable to find any examples for logit regression for non-binary data using Google. I was wondering if there was some reason (e.g. distributional assumptions about the errors) that make it a bad idea to use e.g. glmer with a logit link with continuous data. – Rob Hall Aug 11 '13 at 9:15
• Logit modelling for responses that are continuous proportions has been around for some time, but seemingly is not well known. See e.g. Baum in stata-journal.com/sjpdf.html?articlenum=st0147 Nevertheless it's not your situation. I can't comment on R implementations. – Nick Cox Aug 11 '13 at 9:40
• Thank you Nick for this interesting link - that clears up a few things for me. Sadly, it seems I cannot yet upvote your response. (In case people have trouble accessing the direct link, the following worked for me: stata-journal.com/article.html?article=st0147 ) – Rob Hall Aug 11 '13 at 9:51
• Logistic growth implies a monotonic rising curve. If the data don't match, you'll get a poor fit or the software won't play, depending on exactly what you are doing. – Nick Cox Aug 11 '13 at 9:59

I wanted to share some of the things I learned since asking this question. nlme seems a reasonably way to model non-linear mixed effects in R. Start with a simple base model:

library(nlme)
data <- groupedData(y ~ t | UID, data=data) ## not strictly necessary
initVals <- getInitial(y ~ SSlogis(t, Asym, xmid, scal), data = data)
baseModel<- nlme(y ~ SSlogis(t, Asym, xmid, scal),
data = data,
fixed = list(Asym ~ 1, xmid ~ 1, scal ~ 1),
random = Asym + xmid + scal ~ 1|UID,
start = initVals
)


Then use update to increase model complexity. The start parameter is slightly tricky to work with, it may take some tinkering to figure out the order. Note how the new fixed effect for var1 on Asym follows the regular fixed effect for Asym.

 nestedModel <- update(baseModel, fixed=list(Asym ~ var1, xmid ~ 1, scal ~ 1), start = c(fixef(baseModel)[1], 0, fixef(baseModel)[2], fixef(baseModel)[3]))


lme4 seemed more robust against the outliers in my dataset and seemed to offer more reliable convergence for the more complex models. However, it seems the downside is that the relevant likelihood functions need to be specified manually. The following is the logistic growth model with a fixed effect of var1 (binary) on Asym. You can add fixed effects on xmid and scal in a similar fashion. Note the strange way of specifying the model using a double formula as outcome ~ fixed effects ~ random effects.

library(lme4) ## careful loading nlme and lme4 concurrently
customLogitModel <- function(t, Asym, AsymVar1, xmid, scal) {
(Asym+AsymVar1*var1)/(1+exp((xmid-t)/scal))
}

body(customLogitModel)[[2]],
namevec = c("Asym", "AsymVar1", "xmid", "scal"),
function.arg=customLogitModel
)

## find starting parameters
initVals <- getInitial(y ~ SSlogis(t, Asym, xmid, scal), data = data)

# Fit the model
model <- nlmer(
y ~ customLogitModelGradient(t=t, Asym, AsymVar1, xmid, scal, var1=var) ~
# Random effects with a second ~
(Asym | UID) + (xmid | UID) + (scal | UID),
data = data,
start = c(Asym=initVals[1], AsymVar1=0, xmid=initVals[2], scal=initVals[3])
)

• Thank you Rob for your post, it is actually exactly what I'm trying to do with my data. I don't understand what is var1 in the nestedModel on Asym and how you calculated it ? – user83784 Jul 31 '15 at 9:19
• This is just an example on how to include the effect of some variable on Asym: "The following is the logistic growth model with a fixed effect of var1 (binary) on Asym." For example, you have variable "Treated" that has two values 0 and 1, so substitute "Treated" for "var1". – PA6OTA Mar 16 '17 at 16:58