nlme 错误“组的公式无效”，尽管指定了随机效应

Question

我对此进行了一些搜索，但我发现的邮件列表帖子与未在

nlme

中指定随机效果的人相关，而我已经这样做了。我还拥有 Pinheiro 和 Bates 撰写的《S 和 S-Plus 中的混合效应模型》一书，但无法从书中解决我的问题。

我仍在进行营养数据分析，现在已转向实际数据。这些数据来自一项人口调查，并采用重复测量设计，因为每个受访者都有两次 24 小时内对该营养素的摄入量回忆。

我已经成功地将 lme4 模型拟合到我的数据中，现在我试图找出如果我使用非线性方法会发生什么。我的数据快照如下：

head(Male.Data)
   RespondentID Age SampleWeight  IntakeDay IntakeAmt AgeFactor BoxCoxXY
2        100020  12    0.4952835 Day1Intake 12145.852     9to13 15.61196
7        100419  14    0.3632839 Day1Intake  9591.953    14to18 15.01444
8        100459  11    0.4952835 Day1Intake  7838.713     9to13 14.51458
12       101138  15    1.3258785 Day1Intake 11113.266    14to18 15.38541
14       101214   6    2.1198688 Day1Intake  7150.133      4to8 14.29022
18       101389   5    2.1198688 Day1Intake  5091.528      4to8 13.47928

数据汇总信息为：

str(Male.Data)
'data.frame':   4498 obs. of  7 variables:
$ RespondentID: Factor w/ 4487 levels "100013","100020",..: 2 7 8 12 14 18 19 20 21 22 ...
$ Age         : int  12 14 11 15 6 5 10 2 2 9 ...
$ SampleWeight: num  0.495 0.363 0.495 1.326 2.12 ...
$ IntakeDay   : Factor w/ 2 levels "Day1Intake","Day2Intake": 1 1 1 1 1 1 1 1 1 1 ...
$ IntakeAmt   : num  12146 9592 7839 11113 7150 ...
$ AgeFactor   : Factor w/ 4 levels "1to3","4to8",..: 3 4 3 4 2 2 3 1 1 3 ...
$ BoxCoxXY    : num  15.6 15 14.5 15.4 14.3 ...

使用

lme4

包，我成功地拟合了线性混合效应模型（随机效应来自受试者，

IntakeDay

是与

BoxCoxXY

相关的重复测量因子，它是

IntakeAmt

的变换） :

Male.lme1 <- lmer(BoxCoxXY ~ AgeFactor + IntakeDay + (1|RespondentID),
        data = Male.Data, 
        weights = SampleWeight)

我一直在尝试使用

nlme

包来拟合非线性模型来比较两者，但我无法让我的语法起作用。我最初的问题是，我的数据似乎没有相关的 SelfStart 模型，因此我使用

geeglm

来生成起始值（系数保存到名为

Male.nlme.start

的数据框中）。但现在我收到错误：

Error in getGroups.data.frame(dataMix, eval(parse(text = paste("~1", deparse(groups[[2]]),  : 
Invalid formula for groups

我无法弄清楚我做错了什么，我使用的

nlme

语法是：

Male.nlme1 <- nlme(BoxCoxXY ~ AgeFactor + IntakeDay + RespondentID , data = Male.Data, 
fixed = AgeFactor + IntakeDay ~ 1, 
random = RespondentID ~ 1,
start=c(Male.nlme.start[,"Estimate"]))

我尝试了在整体模型规范中包含和不包含

RespondentID

的分析，这似乎没有影响。

我尝试坚持使用非线性方法的原因是SAS中的原始分析使用了非线性方法。虽然我的残差等从 lme 分析中看起来不错，但我很好奇非线性方法会产生什么影响。

如果有帮助，上次分析尝试的

traceback()

结果（包括

RespondentID

）是：

5: stop("Invalid formula for groups")
4: getGroups.data.frame(dataMix, eval(parse(text = paste("~1", deparse(groups[[2]]), 
   sep = "|"))))
3: getGroups(dataMix, eval(parse(text = paste("~1", deparse(groups[[2]]), 
   sep = "|"))))
2: nlme.formula(BoxCoxXY ~ AgeFactor + IntakeDay, data = Male.Data, 
   fixed = AgeFactor + IntakeDay ~ 1, random = RespondentID ~ 
       1, start = c(Male.nlme.start[, "Estimate"]))
1: nlme(BoxCoxXY ~ AgeFactor + IntakeDay, data = Male.Data, fixed = AgeFactor + 
   IntakeDay ~ 1, random = RespondentID ~ 1, start = c(Male.nlme.start[, 
   "Estimate"])) 

任何人都可以建议我哪里出了问题吗？我开始想知道是否（1）

RespondentID

的因子级别太多而无法在

nlme

中工作，或者（2）只有在我为

RespondentID

提供启动参数时该方法才会起作用，这似乎是无意义的使用我拥有的数据，因为这是我的主题标识符。

更新：回答 Ben，SAS

nlmixed

模型是固定效应的通用对数似然函数：

ll1 <- log(1/sqrt(2*pi*Scale))
ll2 <- as.data.frame(-(BoxCoxXY - Intercept + AgeFactor + IntakeDay + u2)^2)/(2*Scale)+(Lambda.Value-1)*log(IntakeAmt)
ll <- ll1 + ll2
model IntakeAmt ~ general(ll)

地点：

Scale

=

geeglm

和

的色散值

Lambda.Value

= 与早期

boxcox()

输出的最大对数似然相关的 lambda 值，该值用于通过公式

IntakeAmt

 将 

BoxCoxXY

 转换为 

Male.Data$BoxCoxXY <- (Male.Data$IntakeAmt^Lambda.Value-1)/Lambda.Value

SAS代码中的

random

语句是：

random u1 u2 ~ normal([0,0][&vu1,COV_U1U2,&vu2]) subject=RespondentID

因此模型中有两个误差项，它们都被拟合为随机效应。第二个方括号表示按行顺序列出的随机效应方差矩阵的下三角，并使用 SAS 语法中的 SAS 宏变量指定。

我得到的模型摘要是正常的一行概述，显示协变量矩阵 (BX) 加上误差分量，因此这里没有太多帮助。

第二次更新：我意识到我没有删除与女性受试者相关的 RespondentID 级别，因为在按性别将 RespondentID 分解为单独的数据框进行分析之前，我在整个数据框中分解了 RespondentID。在删除 RespondentID 的未使用因子水平后，我重复了

nlme

分析，但出现了相同的错误。

lmer

结果是相同的 - 很高兴知道这一点。 :)

Answer 1

我刚刚偶然发现了这个悬而未决的问题。如果有人仍然对潜在的解决方案感兴趣：

首先，我不是专家。然而，对于一种“非线性”线性回归，我喜欢使用 GAM，例如来自

mgcv

或

gamm4

包。

Wood，S.N.（2017）。广义加性模型：R. chapman 和 Hall/CRC 的介绍。

然后你可以做类似的事情：

library(mgcv)
mygam <- gamm(BoxCoxXY ~ s(AgeFactor) + s(IntakeDay),
             random = list(RespondentID ~ 1),
             data = Male.Data,
             correlation = corAR1(IntakeDay))

请注意，应考虑重复测量（相关观察），例如，在像我在本示例中包含的相关结构中。此外，平滑 (

s()

) 仅适用于连续数据，不适用于随机/离散数据。您可能必须相应地重组数据或包含平滑因子交互（例如：

s(continuous, by = factor)

）。 但是，对于您的明确示例，如果这些数据可能的话，您可能必须调整/重新考虑您的数据结构。

但也许这个提示可以帮助你走上正轨？！