我已经对此进行了一些搜索,但是我发现的邮件列表帖子与未指定随机效果的人相关联,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,SASnlmixed
模型是固定效应的一般对数似然函数:
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()
转换IntakeAmt
为的最大对数似然输出相关联的 lambda 值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 分解为整个数据框之前按性别拆分为单独的数据框以进行分析。nlme
在删除 RespondentID 的未使用因子水平后,我重复了分析,我得到了同样的错误。结果lmer
是一样的——很高兴知道。:)