I am trying to model some data regarding a predator prey interaction experiment (n=26). Predation rate is my response variable and I have 4 explanatory variables: predator density (1,2,3,4 5), predator size, prey density (5,10,15,20,25,30) and prey type (3 categories). I started with several linear models (GLM) and found (as expected) that prey and predator density were non-linearly related to predation rates. If I use a log transformation on these variables I get really nice curves and an adjusted $R^{2}$ of 0.82, but it is not really the right approach for modelling non-linear relationships.

model <-glm(rates ~ log(pred) + log (prey) + type)

Therefore I switched to non-linear least square regression (nls). I have several predator-prey models based on existing ecological literature e.g.:

 ### Holling's type II functional response

model1 <- nls(rates ~ (a * prey)/(1 + b * prey),
start = list(a = 0.27,b = 0.13), trace = TRUE)

### Beddington-DeAngelis functional response

model2 <- nls(rates ~ (a*prey)/(1+ (b * prey) + c * (pred -1 )),
start = list(a=0.22451, b=-0.18938, c=1.06941), trace=TRUE, subset=I1) 

These models work perfectly, but now I want to add prey type as well. In the linear models prey type was the most important variable so I don't want to leave it out. I understand that you can't add categorical variables in nls, so I thought I try a generalized additive model (GAM).

The problem with the gam models is that the smoothers (both spline and loess) don't work on both variables because there are only a very restricted number of values for prey density and predator density. I can manage to get a model with a single variable smoothed using loess. But for two variables it is simply not working. The spline function does not work at all because I have so few values (5) for my variables (see model 4).

model3 <- gam(rates~ lo(pred, span=0.9)+prey)
## this one is actually working but does not include a smoother for prey.

model4 <- gam(rates~ s(pred)+prey)
## this one gives problems: 
A term has fewer unique covariate combinations than specified maximum degrees of freedom

My question is: are there any other possibilities to model data with 2 non-linear related variables in which I can also include a categorical variable. I would prefer to use nls (model2) with for example different intercepts for each category but I'm not sure how to get this sorted, if it is possible at all. The dataset is too small to split it up into the three categories, moreover, one of the categories only contains 5 data points.

  • 1
    $\begingroup$ The gam model starts with knots set at 10, which you do not have enough data for, however you have densities of 5 and 6 groups - you can set your knots 1 fewer and the model should run. Something like model4<-gam(rates~s(pred,k=5)+s(prey,k=4),data=data) $\endgroup$
    – B Williams
    Jul 3, 2013 at 16:20
  • $\begingroup$ @BWilliams Thanks for the advise, the problem with gam was indeed solved. I have also managed to get my head around the nls model. I have created 3 new variables 1 for each category With either the value 1 or 0. I added these in the model as different intercepts. model <- nls(rates ~ fcat1 + dcat2 + e*cat3 +((a * prey)/((1 + b * prey) * (1 + c * (pred-1)))), start = list(a = 0.14, b = 0.009, c=0.66, d=0.8,e=-0.04,f=1.4), trace = TRUE).Not sure if it is legitimate but it works fine $\endgroup$
    – Robbie
    Jul 4, 2013 at 10:45

1 Answer 1


I apologize for the incredibly late response, but I came across the same problem recently. 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))

Although, with 4 different categorical variables, this could get tedious to code. You would have to be very careful with your design matrices.


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