重写混合效果模型公式从R(lme4)到Julia

Question

我想在Julia中使用来自R中lme4库的lmer函数获得相同的结果.请在下面找到使用R的内置mtcars数据集的示例

library(lme4)
data<-mtcars    
data$vs<-as.factor(data$vs)
data$am<-as.factor(data$am)
data$gear<-as.factor(data$gear)
str(data)
model <- lmer(mpg ~ cyl:gear + hp:am + (1|gear:am), data = data)

我已经找到了lmm()来自Julia的MixedModels包的函数,它应该能够运行相同的结果,但是我不知道如何从lmer()函数的第一个参数重写公式lmm().特别是interaction(:)运算符.

我将以简短的例子回答你的问题.

Answer 1

这里 R 和 Julia 模型之间的对应关系似乎并不准确。还有不同数值算法的问题。MixedModels但我尝试使用以下方法重新创建相同的模型：

using RDatasets
using MixedModels

mtcars = dataset("datasets","mtcars")
mtcars[:AM] = PooledDataArray(mtcars[:AM])
mtcars[:Gear] = PooledDataArray(mtcars[:Gear])
mtcars[:GearAM] = PooledDataArray(collect(zip(mtcars[:Gear],mtcars[:AM])))
m = fit!(lmm(MPG ~ 1 + Gear + AM + Cyl + HP + (1|GearAM),mtcars))

手动创建混合效果列很笨拙 - 也许有更好的方法。请注意 R 和 Julia 之间系数命名的差异。两者都有 6 个固定效应系数。

Julia 中的解决方案似乎不同（在我的机器上），但实现了更好的对数似然。随机效应预计会很弱，因为它的变量已经存在于固定效应中（它只考虑了 Gear 和 AM 之间的依赖性），并且只有 32 个数据点。

希望这对您有所帮助，如果您有更好的理解，最好将其添加到另一个答案或评论中。