将混合模型公式从 SAS 转换为 R答案

【问题标题】：Converting mixed model formula from SAS to R将混合模型公式从 SAS 转换为 R
【发布时间】：2011-11-07 01:56:41
【问题描述】：

我想在 R 中使用 nlme 包拟合混合模型，这相当于以下 SAS 代码：

proc mixed data = one;
class var1 var2  year loc rep;
model yld = var1 * var2;
random loc year(loc) rep*year(loc);

编辑：关于什么是实验的解释

var1 和 var2 的相同组合在重复中进行了测试（重复编号为 1:3）。重复（rep）被认为是随机的。这组实验在地点 (loc) 和年份 (year) 上重复。尽管为了方便起见，每个位置和年份内的复制品编号为 1:3，因为它们没有任何名称，但一个位置和一年内的复制品 1 在其他位置和其他年份内没有相关性复制品 1

我尝试了以下代码：

 require(nlme) 
    fm1 <- lme(yld ~ var1*var2, data = one, random = loc + year / loc + rep * year / loc)

我的代码正确吗？

EDITS：基于建议的数据和模型您可以从以下链接下载示例数据文件： https://sites.google.com/site/johndatastuff/mydata1.csv

data$var1 <- as.factor(data$var1)
data$var2 <- as.factor(data$var2)
data$year <- as.factor(data$year)
data$loc <- as.factor(data$loc)
data$rep <- as.factor(data$rep)

following suggestions from the comments below:
fm1 <- lme(yld ~ var1*var2, data = data, random = ~ loc + year / loc + rep * year / loc)

Error in getGroups.data.frame(dataMix, groups) : 
  Invalid formula for groups

基于 SAS 输出的预期

Type 3 tests of fixed effects 
var1*var2         14         238       F value 16.12 Pr >F = < 0.0001

Covariance parameters:
loc = 0, year(loc) = 922161, year*rep(loc) = 2077492, residual = 1109238

我尝试了以下模型，但仍然出现一些错误：

Edits: Just for information I tried the following model
require(lme4)  
 fm1 <- lmer(yld ~ var1*var2 + (1|loc) +  (1|year / loc) + (1|rep : (year / loc)),  
            data = data)  
Error in rep:`:` : NA/NaN argument 
In addition: Warning message: 
In rep:`:` : numerical expression has 270 elements: only the first used

【问题讨论】：

你尝试的时候发生了什么？你能做一个可重复的小例子吗？
您是否尝试在公式中添加~？ fm1
嗯...我没有使用过 lme/lmer 但任何比 lmer(yld ~ var1*var2 +(1|year), data=data) 更高级的东西都失败了。
底部失败是因为“rep : (year”冒号放在一组括号上。这不是 R 所期望的。我不了解 SAS，但似乎你可能想要类似： lmer(yld ~ var1*var2 + (1|loc) + (1|loc/year) + (1|loc/year:rep)), data = data)??
@John 请不要在 SE 网站上发帖。根据经验，纯 R 问题应该继续，统计问题属于这里；否则，请选择一个站点并将您的问题标记为 mod 注意。我正在将其移至 SO，它将与您的其他问题合并以保留 cmets 的历史。（感谢 Aaron 注意到这一点！）

标签： r sas

【解决方案1】：

感谢您提供更详细的信息。我将数据存储在d 中以避免与data 函数和参数混淆；这些命令可以以任何一种方式工作，但这种避免data 通常被认为是好的做法。

请注意，由于var 和var2 之间缺乏平衡，交互很难适应；参考这里的交叉表：

> xtabs(~var1 + var2, data=d)
    var2
var1  1  2  3  4  5
   1 18 18 18 18 18
   2  0 18 18 18 18
   3  0  0 18 18 18
   4  0  0  0 18 18
   5  0  0  0  0 18

通常只适合您使用: 而不是* 的交互（并且没有主效应），但这里最好制作一个单一因素，如下所示：

d$var12 <- factor(paste(d$var1, d$var2, sep=""))

然后用nlme，试试

fm1 <- lme(yld ~ var12, random = ~ 1 | loc/year/rep, data = d)
anova(fm1)

和lme4一起试试

fm1 <- lmer(yld ~ var12 + (1 | loc/year/rep), data=d)
anova(fm1)

还请注意，由于nlme 和lme4 的函数名称有重叠，您只需一次将一个加载到您的 R 会话中；要切换，您需要关闭 R 并重新启动。（还有其他方法，但这是最容易解释的。）

【讨论】：

谢谢；如果你看一下 SAS 模型，我的兴趣主要是交互项而不是主效应，所以即使我也不担心为 var1 和 var2 拟合主效应模型
不平衡是有目的的，1-2之间的预期产品组合与2-1相同，同样1-3与3-1相同。这就是为什么您在 Xtab 的下三角形中找到 0 的原因。
我相信detach(package:lme4) 就足够了（例如，如果lme4 在nlme 之后加载，我们现在想使用nlme::lmList），但我猜你已经知道了.
@chl：实际上，我忘记了如何用另一种方式来做（虽然知道这是可能的），所以谢谢！