Solved – Non-linear mixed effects regression in R

growth-modelmixed modelnonlinear regressionr

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.

Best Answer

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))
}

customLogitModelGradient <- deriv(
    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])
)
Related Question