如何在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)) 

而且這也可以。