如何在R中对lme中的分段混合模型进行编码？

Question

我跟着this example运行分段混合模型使用lmer，它工作得很好。然而，我很难将模型翻译成lme，因为我需要处理异方差，并且lmer没有这种能力。

重现此问题的代码是here。如果您认为有必要回答这个问题，我会在代码中包含有关实验设计的详细信息。

这里是没有断点的模型：

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

这里是我如何与断点运行：

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) 

的问题是，我有很严重的异方差性。据我所知，这意味着我应该使用nlme软件包中的lme。这里是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)) 

这是断点模式，即，好了，破：

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

以下是错误消息：

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

哪有我把这个可爱的断点模型从lmer翻译成lme？谢谢！

Answer 1

看起来像lme不喜欢它，当你在你的公式中使用的变量不在你正在适合你的模型的data.frame中时。一种选择是先建立你的配方，然后将它传递给lme。例如

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

的eval()/substitute()只是在公式中换出bp用的值的变量bp

或者，如果bp总是30，你只想把直接在公式中

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

而且这也可以。