How do I code a piecewise mixed-model in lme in R?

I followed this example for running a piecewise mixed model using lmer, and it works very well. However, I am having trouble translating the model to lme because I need to deal with heteroscedasticity, and lmer doesn’t have that ability.

Code to reproduce the problem is here. I included details about the experimental design in the code if you think it’s necessary to answer the question.

Here is the model without the breakpoint:

linear <- lmer(mass ~ lat + (1 | pop/line), data = df)

And here is how I run it with the breakpoint:

bp = 30
b1 <- function(x, bp) ifelse(x < bp, x, 0)
b2 <- function(x, bp) ifelse(x < bp, 0, x)
breakpoint <- lmer(mass ~ b1(lat, bp) + b2(lat, bp) + (1 | pop/line), data = df)

The problem is that I have pretty severe heteroscedasticity. As far as I understand, that means I should be using lme from the nlme package. Here is the linear model in lme:

ctrl <- lmeControl(opt='optim')
linear2 <- lme(mass ~ lat , random=~1|pop/line, na.action = na.exclude, data=df, control = ctrl, weights=varIdent(form=~1|pop))

And this is the breakpoint model that is, well, breaking:

breakpoint2 <- lme(mass ~ b1(lat, bp) + b2(lat, bp), random=~1|pop/line, na.action = na.exclude, data=df, control = ctrl, weights=varIdent(form=~1|pop))

Here is the error message:

Error in model.frame.default(formula = ~pop + mass + lat + bp + line,  : variable lengths differ (found for 'bp')

How can I translate this lovely breakpoint model from lmer to lme? Thank you!

Solution

Looks like lme doesn't like it when you use variables in your formula that aren't in the data.frame you are fitting your model on. One option would be to build your formula first then pass it to lme. For example

myform <- eval(substitute(mass ~ b1(lat, bp) + b2(lat, bp), list(bp=bp)))
breakpoint2 <- lme(myform, random=~1|pop/line, na.action = na.exclude, data=df, control = ctrl, weights=varIdent(form=~1|pop))

The eval()/substitute() is just to swap out the bp in your formula with the value of the variable bp

Or if bp were always 30, you would just put that directly in the formula

breakpoint2 <- lme(mass ~ b1(lat, 30) + b2(lat, 30), random=~1|pop/line, na.action = na.exclude, data=df, control = ctrl, weights=varIdent(form=~1|pop))

and that would work as well.