是否可以指定没有任何固定效应的 lmer 模型?

机器算法验证 r 混合模式 lme4-nlme
2022-03-19 16:27:27

在 R 中,如何指定没有全局固定效应的 lmer 模型?例如,如果我说类似

lmer(y ~ (1 | group) + (0 + x | group), data = my_df)

拟合模型将是

yij=a+αi+βixij

我如何拟合模型

yij=αi+βixij?

1个回答

正如@Mike Lawrence 提到的,在定义没有固定效果的模型时,显而易见的事情是:

lmer(y ~ -1 + (1|GroupIndicator))

这实际上非常简单;一个定义没有截距或 X 矩阵。这不起作用的基本原因是,正如@maxTC 指出的那样“ lme4 包仅专用于混合模型”。

特别是 lmer() 拟合所做的是通过求解之间的惩罚最小二乘回归来计算轮廓偏差y^y以及球形随机效应u0(方程式(11),参考文献(2))。在计算上,该优化过程利用系统的块结构计算相应系统的 Cholesky 分解(等式(5),参考文献(1))。不设置全局固定效果实际上会以 lmer() 的代码无法应对的方式扭曲该块结构。除其他外,条件期望值u是根据β^的但解决β^问一个从未存在过的矩阵系统的解(矩阵RXX在参考文献(1)中,或LX在参考文献(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));

所以我们要做的是定义一个变量Garbage那只是噪音;和之前一样, lmer() 被指示不使用固定截距,而只使用我们定义的 X 矩阵(在本例中为单列矩阵垃圾)。这个额外的高斯噪声变量将与我们的样本测量误差以及您的随机效应方差不相关。不用说,模型的结构越多,获得不需要但具有统计意义的随机相关性的概率就越小。

所以 lmer() 有一个安慰剂X变量(矩阵)玩,你会得到相关的β为零,您不必以任何方式标准化您的数据(将它们居中,白化它们等)。可能尝试对安慰剂进行几次随机初始化X矩阵也不会受到伤害。“垃圾”的最后一点:使用高斯噪声不是“偶然的”;它在所有方差相等的随机变量中具有最大的熵,因此提供信息增益的机会最小。

显然,这与其说是一种解决方案,不如说是一种计算技巧,但它允许用户有效地指定一个 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!

参考:

  1. DM Bates 和 S. DebRoy 的线性混合模型和惩罚最小二乘法,《多元分析杂志》,第 91 卷第 1 期,2004 年 10 月(链接到预印本
  2. 混合模型的计算方法,Douglas Bates,2012 年 6 月。(链接到源