LMM 的似然比检验得出 P 值为 1？
Likelihood Ratio Test for LMM gives P-Value of 1?
我做了一个实验，在这个实验中，人们必须对个人或非个人的道德困境给出答案。我现在想看看困境的类型和参与者给出的答案（是或否）之间是否存在影响他们反应时间的相互作用。
为此，我使用 lme4 包的 lmer() 函数计算了一个线性混合模型。
我的数据如下所示：
subject condition gender.b age logRT answer dilemma pers_force
1 105 a_MJ1 1 27 5.572154 1 1 1
2 107 b_MJ3 1 35 5.023881 1 1 1
3 111 a_MJ1 1 21 5.710427 1 1 1
4 113 c_COA 0 31 4.990433 1 1 1
5 115 b_MJ3 1 23 5.926926 1 1 1
6 119 b_MJ3 1 28 5.278115 1 1 1
我的函数如下所示：
lmm <- lmer(logRT ~ pers_force * answer + (1|subject) + (1|dilemma),
data = dfb.3, REML = FALSE, control = lmerControl(optimizer="Nelder_Mead"))
以主题和困境作为随机因素。这是输出：
Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: logRT ~ pers_force * answer + (1 | subject) + (1 | dilemma)
Data: dfb.3
Control: lmerControl(optimizer = "Nelder_Mead")
AIC BIC logLik deviance df.resid
-13637.3 -13606.7 6825.6 -13651.3 578
Scaled residuals:
Min 1Q Median 3Q Max
-3.921e-07 -2.091e-07 2.614e-08 2.352e-07 6.273e-07
Random effects:
Groups Name Variance Std.Dev.
subject (Intercept) 3.804e-02 1.950e-01
dilemma (Intercept) 0.000e+00 0.000e+00
Residual 1.155e-15 3.398e-08
Number of obs: 585, groups: subject, 148; contrasts, 4
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.469e+00 1.440e-02 379.9
pers_force1 -1.124e-14 5.117e-09 0.0
answer -1.095e-15 4.678e-09 0.0
pers_force1:answer -3.931e-15 6.540e-09 0.0
Correlation of Fixed Effects:
(Intr) prs_f1 answer
pers_force1 0.000
answer 0.000 0.447
prs_frc1:aw 0.000 -0.833 -0.595
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see ?isSingular
然后我使用简化模型进行了似然比检验以获得 p 值：
lmm_null <- lmer(logRT ~ pers_force + answer + (1|subject) + (1|dilemma),
data = dfb.3, REML = FALSE,
control = lmerControl(optimizer="Nelder_Mead"))
anova(lmm,lmm_null)
对于这两个模型，我都收到警告“边界（奇异）拟合：参见 ?isSingular”，但如果我删除一个随机效应以使结构更简单，则会收到模型无法收敛的警告（这有点奇怪），所以我忽略了它。
但是，LRT 输出如下所示：
Data: dfb.3
Models:
lmm_null: logRT ~ pers_force + answer + (1 | subject) + (1 | dilemma)
lmm: logRT ~ pers_force * answer + (1 | subject) + (1 | dilemma)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
lmm_null 6 -13639 -13613 6825.6 -13651
lmm 7 -13637 -13607 6825.6 -13651 0 1 1
可以看到，卡方值为0，p-Value正好为1，看起来很奇怪。我想一定是这里出了什么问题，但我不知道是什么。
你说
logRT is the average logarithmized reaction time across all those 4 dilemmas.
如果我对此的解释是正确的——即每个受试者在观察他们的所有时间里都有相同的反应——那么这就是你的问题的近端原因。 （我知道我以前见过这个确切的问题，但我不知道在哪里 — 这里？r-sig-mixed-models@r-project.org？）
模拟数据
library(lme4)
set.seed(101)
dd1 <- expand.grid(subject=factor(100:150), contrasts=factor(1:4))
dd1$answer <- rbinom(nrow(dd1),size=1,prob=0.5)
dd1$logRT <- simulate(~answer + contrasts + (1|subject),
family=gaussian,
newparams=list(beta=c(0,1,1,-1,2),theta=1,sigma=1),
newdata=dd1)[[1]]
常规合身
这很好，给出了接近真实参数的答案：
m1 <- lmer(logRT~answer + contrasts + (1|subject), data=dd1)Linear mixed model fit by REML ['lmerMod']
## Random effects:
## Groups Name Std.Dev.
## subject (Intercept) 1.0502
## Residual 0.9839
## Number of obs: 204, groups: subject, 51
## Fixed Effects:
## (Intercept) answer contrasts2 contrasts3 contrasts4
## -0.04452 0.85333 1.16785 -1.07847 1.99243
现在按主题平均回复
我们收到大量警告消息，以及您所看到的相同病状（残差和截距以外的所有参数估计实际上为零）。这是因为 lmer 正在尝试从 within-subject 变化中估计残差，我们已经摆脱了它！
我不知道你为什么要取平均。如果这是不可避免的，并且您的设计是此处显示的随机块类型（每个受试者看到所有四个dilemmas/contrasts），那么您无法估计困境效应。
dd2 <- transform(dd1, logRT=ave(logRT,subject))
m2 <- update(m1, data=dd2)
## Random effects:
## Groups Name Std.Dev.
## subject (Intercept) 6.077e-01
## Residual 1.624e-05
## Number of obs: 204, groups: subject, 51
## Fixed Effects:
## (Intercept) answer contrasts2 contrasts3 contrasts4
## 9.235e-01 1.031e-10 -1.213e-11 -1.672e-15 -1.011e-11
将困境视为随机效应不会如您所愿（考虑到它们的呈现方式因人而异）。困境中的受试者间变异性将被归并到它所属的剩余变异性中——我建议将其视为固定效应。
我做了一个实验，在这个实验中，人们必须对个人或非个人的道德困境给出答案。我现在想看看困境的类型和参与者给出的答案（是或否）之间是否存在影响他们反应时间的相互作用。
为此，我使用 lme4 包的 lmer() 函数计算了一个线性混合模型。
我的数据如下所示：
subject condition gender.b age logRT answer dilemma pers_force
1 105 a_MJ1 1 27 5.572154 1 1 1
2 107 b_MJ3 1 35 5.023881 1 1 1
3 111 a_MJ1 1 21 5.710427 1 1 1
4 113 c_COA 0 31 4.990433 1 1 1
5 115 b_MJ3 1 23 5.926926 1 1 1
6 119 b_MJ3 1 28 5.278115 1 1 1
我的函数如下所示：
lmm <- lmer(logRT ~ pers_force * answer + (1|subject) + (1|dilemma),
data = dfb.3, REML = FALSE, control = lmerControl(optimizer="Nelder_Mead"))
以主题和困境作为随机因素。这是输出：
Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: logRT ~ pers_force * answer + (1 | subject) + (1 | dilemma)
Data: dfb.3
Control: lmerControl(optimizer = "Nelder_Mead")
AIC BIC logLik deviance df.resid
-13637.3 -13606.7 6825.6 -13651.3 578
Scaled residuals:
Min 1Q Median 3Q Max
-3.921e-07 -2.091e-07 2.614e-08 2.352e-07 6.273e-07
Random effects:
Groups Name Variance Std.Dev.
subject (Intercept) 3.804e-02 1.950e-01
dilemma (Intercept) 0.000e+00 0.000e+00
Residual 1.155e-15 3.398e-08
Number of obs: 585, groups: subject, 148; contrasts, 4
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.469e+00 1.440e-02 379.9
pers_force1 -1.124e-14 5.117e-09 0.0
answer -1.095e-15 4.678e-09 0.0
pers_force1:answer -3.931e-15 6.540e-09 0.0
Correlation of Fixed Effects:
(Intr) prs_f1 answer
pers_force1 0.000
answer 0.000 0.447
prs_frc1:aw 0.000 -0.833 -0.595
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see ?isSingular
然后我使用简化模型进行了似然比检验以获得 p 值：
lmm_null <- lmer(logRT ~ pers_force + answer + (1|subject) + (1|dilemma),
data = dfb.3, REML = FALSE,
control = lmerControl(optimizer="Nelder_Mead"))
anova(lmm,lmm_null)
对于这两个模型，我都收到警告“边界（奇异）拟合：参见 ?isSingular”，但如果我删除一个随机效应以使结构更简单，则会收到模型无法收敛的警告（这有点奇怪），所以我忽略了它。 但是，LRT 输出如下所示：
Data: dfb.3
Models:
lmm_null: logRT ~ pers_force + answer + (1 | subject) + (1 | dilemma)
lmm: logRT ~ pers_force * answer + (1 | subject) + (1 | dilemma)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
lmm_null 6 -13639 -13613 6825.6 -13651
lmm 7 -13637 -13607 6825.6 -13651 0 1 1
可以看到，卡方值为0，p-Value正好为1，看起来很奇怪。我想一定是这里出了什么问题，但我不知道是什么。
你说
logRT is the average logarithmized reaction time across all those 4 dilemmas.
如果我对此的解释是正确的——即每个受试者在观察他们的所有时间里都有相同的反应——那么这就是你的问题的近端原因。 （我知道我以前见过这个确切的问题，但我不知道在哪里 — 这里？r-sig-mixed-models@r-project.org？）
模拟数据
library(lme4)
set.seed(101)
dd1 <- expand.grid(subject=factor(100:150), contrasts=factor(1:4))
dd1$answer <- rbinom(nrow(dd1),size=1,prob=0.5)
dd1$logRT <- simulate(~answer + contrasts + (1|subject),
family=gaussian,
newparams=list(beta=c(0,1,1,-1,2),theta=1,sigma=1),
newdata=dd1)[[1]]
常规合身
这很好，给出了接近真实参数的答案：
m1 <- lmer(logRT~answer + contrasts + (1|subject), data=dd1)Linear mixed model fit by REML ['lmerMod']
## Random effects:
## Groups Name Std.Dev.
## subject (Intercept) 1.0502
## Residual 0.9839
## Number of obs: 204, groups: subject, 51
## Fixed Effects:
## (Intercept) answer contrasts2 contrasts3 contrasts4
## -0.04452 0.85333 1.16785 -1.07847 1.99243
现在按主题平均回复
我们收到大量警告消息，以及您所看到的相同病状（残差和截距以外的所有参数估计实际上为零）。这是因为 lmer 正在尝试从 within-subject 变化中估计残差，我们已经摆脱了它！
我不知道你为什么要取平均。如果这是不可避免的，并且您的设计是此处显示的随机块类型（每个受试者看到所有四个dilemmas/contrasts），那么您无法估计困境效应。
dd2 <- transform(dd1, logRT=ave(logRT,subject))
m2 <- update(m1, data=dd2)
## Random effects:
## Groups Name Std.Dev.
## subject (Intercept) 6.077e-01
## Residual 1.624e-05
## Number of obs: 204, groups: subject, 51
## Fixed Effects:
## (Intercept) answer contrasts2 contrasts3 contrasts4
## 9.235e-01 1.031e-10 -1.213e-11 -1.672e-15 -1.011e-11
将困境视为随机效应不会如您所愿（考虑到它们的呈现方式因人而异）。困境中的受试者间变异性将被归并到它所属的剩余变异性中——我建议将其视为固定效应。