我正在使用具有逻辑链接功能的混合效果模型（在 R 中使用 lme4 版本 1.1-7）。但是，我注意到固定效应的显着性估计值会根据数据集中行的顺序而变化。
也就是说，如果我在数据集上运行模型，我会得到对我的固定效应的一定估计，并且它具有一定的 p 值。我再次运行模型，得到相同的估计值和 p 值。现在，我打乱了行的顺序（数据没有混合，只是行的顺序不同）。第三次运行模型，p 值非常不同。
对于我拥有的数据，固定效应的估计 p 值可以在 p=0.001 和 p=0.08 之间。显然，考虑到传统的显着性水平，这些是至关重要的差异。
我了解这些估计值只是估计值，由于多种原因，值之间会存在差异。但是，我的数据差异幅度对我来说似乎很大，我不希望我的数据帧的顺序会产生这种影响（当同事运行相同的模型但得到不同的结果时，我们偶然发现了这个问题。它原来他们已经订购了他们的数据框。）。
这是我的脚本的输出：（X 和 Y 是对比编码和居中的二进制变量，Group 和 SubGroup 是分类变量）
> # Fit model
> m1 = glmer(X ~Y+(1+Y|Group)+(1+Y|SubGroup),family=binomial(link='logit'),data=d)
> # Shuffle order of rows
> d = d[sample(1:nrow(d)),]
> # Fit model again
> m2 = glmer(X ~Y+(1+Y|Group)+(1+Y|SubGroup),family=binomial(link='logit'),data=d)
> summary(m1)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation)
['glmerMod']
Family: binomial ( logit )
Formula: X ~ Y + (1 + Y | Group) + (1 + Y | SubGroup)
Data: d
AIC BIC logLik deviance df.resid
200692.0 200773.2 -100338.0 200676.0 189910
Scaled residuals:
Min 1Q Median 3Q Max
-1.1368 -0.5852 -0.4873 -0.1599 6.2540
Random effects:
Groups Name Variance Std.Dev. Corr
SubGroup (Intercept) 0.2939 0.5421
Y1 0.1847 0.4298 -0.79
Group (Intercept) 0.2829 0.5319
Y1 0.4640 0.6812 -0.07
Number of obs: 189918, groups: SubGroup, 15; Group, 12
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.0886 0.1325 -8.214 <2e-16 ***
Y1 0.3772 0.2123 1.777 0.0756 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
Y1 0.112
>
> # -----------------
> summary(m2)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation)
['glmerMod']
Family: binomial ( logit )
Formula: X ~ Y + (1 + Y | Group) + (1 + Y | SubGroup)
Data: d
AIC BIC logLik deviance df.resid
200692.0 200773.2 -100338.0 200676.0 189910
Scaled residuals:
Min 1Q Median 3Q Max
-1.1368 -0.5852 -0.4873 -0.1599 6.2540
Random effects:
Groups Name Variance Std.Dev. Corr
SubGroup (Intercept) 0.2939 0.5422
Y1 0.1846 0.4296 -0.79
Group (Intercept) 0.2829 0.5318
Y1 0.4641 0.6813 -0.07
Number of obs: 189918, groups: SubGroup, 15; Group, 12
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.0886 0.1166 -9.334 < 2e-16 ***
Y1 0.3773 0.1130 3.339 0.000841 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
Y1 0.074
由于隐私原因，我担心无法附加数据。
两种模型都收敛。差异似乎在于标准误差，而系数估计的差异较小。模型拟合（AIC 等）是相同的，因此可能存在多个最优收敛，并且数据的顺序将优化器推入不同的收敛。但是，每次对数据帧进行洗牌时，我都会得到略微不同的估计（不仅仅是两个或三个唯一的估计）。在一种情况下（上面未显示），模型不会仅仅因为行的洗牌而收敛。
我怀疑问题出在我的特定数据的结构上。它相当大（近 200,000 个案例），并且具有嵌套的随机效应。我已经尝试将数据居中，使用对比编码并根据先前的拟合将起始值提供给 lmer。这似乎有所帮助，但我仍然得到相当大的 p 值差异。我也尝试使用不同的方法来计算 p 值，但我遇到了同样的问题。
下面，我尝试用合成数据复制这个问题。这里的差异没有我的真实数据那么大，但它给出了问题的概念。
library(lme4)
set.seed(999)
# make a somewhat complex data frame
x = c(rnorm(10000),rnorm(10000,0.1))
x = sample(x)
y = jitter(x,amount=10)
a = rep(1:20,length.out=length(x))
y[a==1] = jitter(y[a==1],amount=3)
y[a==2] = jitter(x[a==2],amount=1)
y[a>3 & a<6] = rnorm(sum(a>3 & a<6))
# convert to binary variables
y = y >0
x = x >0
# make a data frame
d = data.frame(x1=x,y1=y,a1=a)
# run model
m1 = glmer(x1~y1+(1+y1|a1),data=d,family=binomial(link='logit'))
# shuffle order of rows
d = d[sample(nrow(d)),]
# run model again
m2 = glmer(x1~y1+(1+y1|a1),data=d,family=binomial(link='logit'))
# show output
summary(m1)
summary(m2)
一种解决方案是使用不同的行顺序多次运行模型，并报告 p 值的范围。然而，这似乎不优雅，并且可能相当混乱。
该问题不会影响模型比较估计（使用方差分析），因为这些是基于模型拟合的差异。固定效应系数估计也相当稳健。因此，我可以只报告效应大小、置信区间和模型与空模型比较的 p 值，而不是主模型内的 p 值。
无论如何，其他人有这个问题吗？关于如何进行的任何建议？