I have a question regarding the use of a grouping variable in a non-linear model. Since the nls() function does not allow for factor variables, I have been struggling to figure out if one can test the effect of a factor on the model fit. I have included an example below where I want to fit a "seasonalized von Bertalanffy" growth model to different growth treatments (most commonly applied to fish growth). I would like to test the effect of the lake where the fish grew as well as the food given (just an artificial example). I am familiar with a workaround to this problem - applying an F-test comparing models fit to pooled data vs. separate fits as outlined by Chen et al. (1992) (ARSS - "Analysis of residual sum of squares"). In other words, for the example below, does the fitting of two models significantly reduce sum of the squared residuals (in this example, yes):

enter image description here

I imagine there is a simpler way to do this in R using nlme(), but I am running into problems. First of all, by using a grouping variable, the degrees of freedom is higher than I obtain with my fitting of separate models. Second, I am unable to nest grouping variables - I don't see where my problem is. Any help using nlme or other methods is greatly appreciated. Below is code for my artificial example:

###seasonalized von Bertalanffy growth model
soVBGF <- function(S.inf, k, age, age.0, age.s, c){
    S.inf * (1-exp(-k*((age-age.0)+(c*sin(2*pi*(age-age.s))/2*pi)-(c*sin(2*pi*(age.0-age.s))/2*pi))))

###Make artificial data
food <- c("corn", "corn", "wheat", "wheat")
lake <- c("king", "queen", "king", "queen")

#cornking, cornqueen, wheatking, wheatqueen
S.inf <- c(140, 140, 130, 130)
k <- c(0.5, 0.6, 0.8, 0.9)
age.0 <- c(-0.1, -0.05, -0.12, -0.052)
age.s <- c(0.5, 0.5, 0.5, 0.5)
cs <- c(0.05, 0.1, 0.05, 0.1)

PARS <- data.frame(food=food, lake=lake, S.inf=S.inf, k=k, age.0=age.0, age.s=age.s, c=cs)

#make data
db <- c()
PCH <- NaN*seq(4)
COL <- NaN*seq(4)
for(i in seq(4)){
    age <- runif(min=0.2, max=5, 100)
    age <- age[order(age)]
    size <- soVBGF(PARS$S.inf[i], PARS$k[i], age, PARS$age.0[i], PARS$age.s[i], PARS$c[i]) + rnorm(length(age), sd=3)
	PCH[i] <- c(1,2)[which(levels(PARS$food) == PARS$food[i])]
	COL[i] <- c(2,3)[which(levels(PARS$lake) == PARS$lake[i])]
	db <- rbind(db, data.frame(age=age, size=size, food=PARS$food[i], lake=PARS$lake[i], pch=PCH[i], col=COL[i]))

#visualize data
plot(db$size ~ db$age, col=db$col, pch=db$pch)
legend("bottomright", legend=paste(PARS$food, PARS$lake), col=COL, pch=PCH)

###fit growth model

starting.values <- c(S.inf=140, k=0.5, c=0.1, age.0=0, age.s=0)

#fit to pooled data ("small model")
fit0 <- nls(size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 

#fit to each lake separatly ("large model")
fit.king <- nls(size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 

fit.queen <- nls(size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 

#analysis of residual sum of squares (F-test)
resid.small <- resid(fit0)
resid.big <- c(resid(fit.king),resid(fit.queen))
df.small <- summary(fit0)$df
df.big <- summary(fit.king)$df+summary(fit.queen)$df

F.value <- ((sum(resid.small^2)-sum(resid.big^2))/(df.big[1]-df.small[1])) / (sum(resid.big^2)/(df.big[2]))
P.value <- pf(F.value , (df.big[1]-df.small[1]), df.big[2], lower.tail = FALSE)
F.value; P.value

###plot models
plot(db$size ~ db$age, col=db$col, pch=db$pch)
legend("bottomright", legend=paste(PARS$food, PARS$lake), col=COL, pch=PCH)
legend("topleft", legend=c("soVGBF pooled", "soVGBF king", "soVGBF queen"), col=c(1,2,3), lwd=2)

#plot "small" model (pooled data)
tmp <- data.frame(age=seq(min(db$age), max(db$age),,100))
pred <- predict(fit0, tmp)
lines(tmp$age, pred, col=1, lwd=2)

#plot "large" model (seperate fits)
tmp <- data.frame(age=seq(min(db$age), max(db$age),,100), lake="king")
pred <- predict(fit.king, tmp)
lines(tmp$age, pred, col=2, lwd=2)
tmp <- data.frame(age=seq(min(db$age), max(db$age),,100), lake="queen")
pred <- predict(fit.queen, tmp)
lines(tmp$age, pred, col=3, lwd=2)

###Can this be done in one step using a grouping variable?
#with "lake" as grouping variable
starting.values <- c(S.inf=140, k=0.5, c=0.1, age.0=0, age.s=0)
fit1 <- nlme(model = size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 
  fixed = S.inf + k + c + age.0 + age.s ~ 1,
  group = ~ lake,

#similar residuals to the seperatly fitted models

#but different degrees of freedom? (10 vs. 21?)
AIC(fit1, fit0)

###I would also like to nest my grouping factors. This doesn't work...
#with "lake" and "food" as grouping variables
starting.values <- c(S.inf=140, k=0.5, c=0.1, age.0=0, age.s=0)
fit2 <- nlme(model = size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 
  fixed = S.inf + k + c + age.0 + age.s ~ 1,
  group = ~ lake/food,

Reference: Chen, Y., Jackson, D.A. and Harvey, H.H., 1992. A comparison of von Bertalanffy and polynomial functions in modelling fish growth data. 49, 6: 1228-1235.


2 Answers 2


You could stratify by the values of the categorical predictor and compare fits. For example suppose you have continuous predictors $X_{1}, ..., X_{p}$ and dependent variable $Y$. I believe nls() gives the maximum likelihood estimate of $f$ such that

$$ Y = f(X_1, ..., X_p) + \varepsilon $$

where $\varepsilon \sim N(0,\sigma^2)$ and $f$ is parameterized in some non-linear way (see the nls helpfile). Suppose you have a categorical predictor $B$ with $m$ levels and stratify the data based on the values of $B$ and fit the model within each strata. Since these are disjoint subsets of the data, the log-likelihood for the stratified model, $L_1$ is the sum of the likelihood within each strata, which can be extracted from an nls model using the logLik function (you can also get the log-likelihood from the unstratified model, $L_0$, using logLik).

The unstratified model is clearly a submodel of the stratified model, so the likelihood ratio test is appropriate to see whether the larger model is worth the added complexity - the test statistic is

$$ \lambda = 2(L_1-L_0)$$

If the categorical predictor truly has no effect, $\lambda$ has a $\chi^2$ distribution with degrees of freedom equal to $mp - p = p(m-1)$, where $p$ is the number of parameters underlying your non-linear regression function. You can use pchisq() to calculate $\chi^2$ p-values.

  • $\begingroup$ Are you suggesting to fit m separate models, sum the log likelihood from each L1= SUM(LL_i, i from 1 to m) and then proceed with the likelihood? Also, is L0 a model with the categorical predictor in question included (with m-1 dummy variables for example)? $\endgroup$
    – B_Miner
    Commented Apr 18, 2012 at 1:48
  • $\begingroup$ Yes, I am suggesting that. $L_0$ is the maximized likelihood when you've left $B$ out of the model entirely (i.e. you're fitting the regression function for the entire data set, not stratifying by values of $B$). $\endgroup$
    – Macro
    Commented Apr 18, 2012 at 2:21
  • $\begingroup$ Thanks for your suggestion Macro. This seems to be in the direction of what I have already done - although you suggest comparison of likelihood rather than the F-test. In my example, the F-test also compares the single fit residuals to the summation of residuals from several fits applied to each categorical predictor level. I guess I was wondering if one can do this within a mixed model in one step rather than fitting several models. Also, would such a strategy allow for nested factor testing? $\endgroup$ Commented Apr 18, 2012 at 5:20
  • $\begingroup$ I don't think you'll be able to get around fitting several models in order to compare models. Also, yes, the likelihood ratio test can be used to test for nested factors. $\endgroup$
    – Macro
    Commented Apr 30, 2012 at 2:25

I found that it is possible to code categorical variables with nls(), simply by multiplying true/false vectors into your equation. Example:

# null model (no difference between groups; all have the same coefficients)
nls.null <- nls(formula = percent_on_cells ~ vmax*(Time/(Time+km)),
            data = mehg,
            start = list(vmax = 0.6, km = 10))

# alternative model (each group has different coefficients)
nls.alt <- nls(formula = percent_on_cells ~ 
            + as.numeric(DOC==1)*(vmax2)*(Time/(Time+(km2)))
            + as.numeric(DOC==10)*(vmax3)*(Time/(Time+(km3)))
            + as.numeric(DOC==100)*(vmax4)*(Time/(Time+(km4))),
            data = mehg, 
            start = list(vmax1=0.63, km1=3.6, 
                         vmax2=0.64, km2=3.6, 
                         vmax3=0.50, km3=3.2,
                         vmax4= 0.40, km4=9.7))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.