r - LMM 的似然比检验给出的 P 值为 1?
问题描述
我做了一个实验,在这个实验中,人们必须对个人或非个人的道德困境给出答案。我现在想看看困境的类型和参与者给出的答案(是或否)之间是否存在相互作用,从而影响他们的反应时间。为此,我使用lmer()
lme4-package 的 -function 计算了一个线性混合模型。我的数据如下所示:
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
可以看到,Chi-Square 值是 0,而 p-Value 正好是 1,这看起来很奇怪。我想这里一定出了什么问题,但我不知道是什么。
解决方案
你说
logRT 是所有这 4 个困境的平均对数反应时间。
如果我对此的解释是正确的——即,每个受试者在被观察的所有时间里都有相同的反应——那么这就是你的问题的近因。(我知道我以前见过这个确切的问题,但我不知道在哪里——在这里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
试图从主体内的变异中估计残差,而我们已经摆脱了它!
我不知道你为什么要进行平均。如果这是不可避免的,并且您的设计是此处显示的随机区组类型(每个受试者都看到所有四个困境/对比),那么您无法估计困境效应。
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
将困境视为随机效应不会达到您想要的效果(允许个体与个体之间的差异)。困境中的主体间变异性将被归入它所属的剩余变异性中——我建议将其视为固定效应。
推荐阅读
- ios - 如何将参数中的 ID 发送到我的服务器?
- php - CI_DB_mysqli_result 类的对象无法转换为字符串?
- c - 标准输入和标准输出实际上是同一个文件吗?
- cordova - 日期选择器在更改事件上打开键盘
- azure-managed-identity - 如何在 VM/VMSS 上为特定用户分配的身份获取 Azure MSI 访问令牌?
- android - android.support.v4.media.MediaBrowserCompat$CustomActionCallback AndroidX
- c# - 当页面超过一页时,Rotativa 会破坏 HTML
- hadoop - 将多个 Hive 表合并为 Hadoop 中的单个表
- wordpress - Wordpress WP_User_Query + Relevanssi
- javascript - 文件阅读器 Angular 2