在 R 中,如何指定没有全局固定效应的 lmer 模型?例如,如果我说类似
lmer(y ~ (1 | group) + (0 + x | group), data = my_df)
拟合模型将是
我如何拟合模型
?
在 R 中,如何指定没有全局固定效应的 lmer 模型?例如,如果我说类似
lmer(y ~ (1 | group) + (0 + x | group), data = my_df)
拟合模型将是
我如何拟合模型
?
正如@Mike Lawrence 提到的,在定义没有固定效果的模型时,显而易见的事情是:
lmer(y ~ -1 + (1|GroupIndicator))
这实际上非常简单;一个定义没有截距或 X 矩阵。这不起作用的基本原因是,正如@maxTC 指出的那样“ lme4 包仅专用于混合模型”。
特别是 lmer() 拟合所做的是通过求解之间的惩罚最小二乘回归来计算轮廓偏差和以及球形随机效应 和(方程式(11),参考文献(2))。在计算上,该优化过程利用系统的块结构计算相应系统的 Cholesky 分解(等式(5),参考文献(1))。不设置全局固定效果实际上会以 lmer() 的代码无法应对的方式扭曲该块结构。除其他外,条件期望值是根据的但解决问一个从未存在过的矩阵系统的解(矩阵在参考文献(1)中,或在参考文献(2))。因此,您会收到如下错误:
Error in mer_finalize(ans) :
Cholmod error 'invalid xtype' at file:../Cholesky/cholmod_solve.c, line 970
因为毕竟一开始就没有什么可以解决的。
假设您不想重新编写 lmer() 分析的偏差成本函数,最简单的解决方案是基于 CS-101 公理:垃圾输入,垃圾输出。
N = length(y); Garbage <- rnorm(N);
lmer(y ~ -1 + Garbage + (1|GroupIndicator));
所以我们要做的是定义一个变量那只是噪音;和之前一样, lmer() 被指示不使用固定截距,而只使用我们定义的 X 矩阵(在本例中为单列矩阵垃圾)。这个额外的高斯噪声变量将与我们的样本测量误差以及您的随机效应方差不相关。不用说,模型的结构越多,获得不需要但具有统计意义的随机相关性的概率就越小。
所以 lmer() 有一个安慰剂变量(矩阵)玩,你会得到相关的为零,您不必以任何方式标准化您的数据(将它们居中,白化它们等)。可能尝试对安慰剂进行几次随机初始化矩阵也不会受到伤害。“垃圾”的最后一点:使用高斯噪声不是“偶然的”;它在所有方差相等的随机变量中具有最大的熵,因此提供信息增益的机会最小。
显然,这与其说是一种解决方案,不如说是一种计算技巧,但它允许用户有效地指定一个 lmer 模型而没有全局固定效应。很抱歉希望围绕这两个参考。总的来说,我认为 Ref.(1) 是任何人了解 lmer() 正在做什么的最佳选择,但 Ref.(2) 更接近实际代码的精神。
这里有一些代码展示了上面的想法:
library(lme4)
N= 500; #Number of Samples
nlevA = 25; #Number of levels in the random effect
set.seed(0) #Set the seed
e = rnorm(N); e = 1*(e - mean(e) )/sd(e); #Some errors
GroupIndicator = sample(nlevA, N, replace=T) #Random Nvel Classes
Q = lmer( rnorm(N) ~ (1| GroupIndicator )); #Dummy regression to get the matrix Zt easily
Z = t(Q@Zt); #Z matrix
RA <- rnorm(nlevA ) #Random Normal Matrix
gammas =c(3*RA/sd(RA)) #Colour this a bit
y = as.vector( Z %*% gammas + e ) #Our measurements are the sum of measurement error (e) and Group specific variance
lmer_native <- lmer(y ~ -1 +(1| GroupIndicator)) #No luck here.
Garbage <- rnorm(N) #Prepare the garbage
lmer_fooled <- lmer(y ~ -1 + Garbage+(1| GroupIndicator)) #OK...
summary(lmer_fooled) #Hey, it sort of works!
参考: