基于线性混合效应模型的心理实验数据建模

/基于线性混合效应模型的心理实验数据建模【摘要】线性混合效应模型(LMEM)作为一种数据分析方法在心理科学中越来越受欢迎。在心理实验的应用中，线性混合效应模型仍囿于正确反应时，而错误反应时也反映了被试的重要操作信息。基于此，借鉴生存分析对截尾数据的处理思路，文中提出联合错误反应时的LMEM(JILMEM)。在实例数据中，JILMEM具有更佳的模型拟合。在模拟研究中，相对LMEM和F1/F2检验方法，JILMEM对不同模拟情景特点和处理效应识别更有优势。因此，JILMEM能更有效识别实验设计特点和更佳性能。 研究要点 1.基于生存分析提出了可联合分析错误反应时的线性混合效应模型。 2.基于实例与模拟数据系统评估不同模型与方法的性能。 1引言 心理学实验常通过样本结果推断出总体的概括性结论(Barr，2017)。然而，传统的方差分析无法同时考虑实验设计中被试与刺激的随机特点。例如，常用的F1/F2检验通过计算被试平均反应时或刺激平均反应时后进行方差分析(吴岩，李天虹，2018)，无论哪种做法都存在“随机效应固定化”的问题(Juddetal.，2012)。因此，线性混合效应模型(LinearMixed-EffectsModel，LMEM)因具有灵活处理多个随机变量的优势而受到心理实验研究者关注(Baayenetal.，2008)，逐渐成为心理学和语言学等实验研究的数据分析方法的重要新选择(Meteyard&Davies，2020)。 LMEM最早由Baayen等(2008)引入心理语言学实验数据分析中，随后不同研究者从不同方面展开了研究。在LMEM的随机效应结构的探索上，Judd等人(2012)提出了被试和刺激作为随机变量的理论随机效应结构；同时，Judd等人(2017)对不同实验设计的LMEM的随机效应结构进行了理论探讨。在随机效应结构的最优选择上，Barr等人(2013)通过模拟发现，在标准误和第I类错误率方面，包含实验设计的最大随机效应结构表现最好。然而，Matuschek等人(2017)通过数据驱动来确定随机效应结构。Brauer和Curtin(2018)进一步探讨了在LMEM中整合连续和类别自变量的随机效应结构特点。对LMEM的参数估计，不同研究者分别基于R语言包、JAGS和stan(Brown，2021；Sorensenetal.，2016)介绍相应估计方法。对于LMEM的使用指导，一些研究提供了统计检验力和效应量计算的教程以及实践应用指导(Brysbaert&Stevens，2018；DeBruine&Barr，2021；Meteyardetal.，2020；马拯等，2022)。然而，现有的LMEM仍然囿于正确反应时，常采用替换或删除的方式处理错误反应时的数据。错误反应时也反映了被试在刺激上的重要操作信息(Milleretal.，2018)，此时现有的LMEM并不适用。 对于错误反应时信息，Miller等(2018)认为错误反应时是正确反应时的截尾点。在生存分析中，截尾是指由于时间的限制而没有观测到某些事件发生的时间，常通过生存函数描述截尾数据发生的可能性(Miller，1998)。因此，现有的LMEM通过平均正确反应时替换或者删除错误反应时的方法可能并不恰当，开发相应的LMEM是有必要的。为了解决上述问题，发挥LMEM的优势，本研究借鉴生存分析对截尾数据的处理思路，提出一个联合错误反应时的LMEM，即JILMEM(TheLMEMofJointIncorrectReactionTime)，并综合比较不同数据分析方法，以期为心理学实验数据分析提供有价值的指导建议。 2联合错误反应时的LMEM 2.1模型构建 在介绍JILMEM前，先简单介绍LMEM。在不同实验设计中，因被试与刺激安排不同，随机效应结构也是不相同的(Juddetal.，2017)。本研究聚焦最常用的单因素两水平重复测量实验设计(DeBruineetal.，2021)，对LMEM展开研究。设定Tijk为被试i(i=1，…，N)在条件k(k=1，…，K)对刺激项目j(j=1，…，M)的反应时，则LMEM可以表示为 其中，β0和β1为固定效应，前者为总体平均反应时，后者为不同实验条件Xk斜率或处理效应，Xk为自变量不同水平，为比较不同水平差异和估计模型截距，需要编码为虚拟变量，常见的两水平编码为X1=0.5和X2=-0.5(DeBruineetal.，2021；Matuscheketal.，2017)。S0i和S1i是被试水平的随机截距和斜率，假设服从二元正态分布，其均值向量为0，方差分别为τ200和τ211，相关系数为ρs(式2)；I0i是刺激项目水平的随机截距，因刺激项目嵌套于不同实验条件，并不存在刺激项目水平的随机斜率(Juddetal.，2012)，I0i假设服从正态分布N(0，ω200)；ξijk为残差反应时，假设服从均值为0，方差为σ2的正态分布。 为了更好地理解JILMEM的构建思路，以及比较LMEM和JILMEM的差异，将从认知过程模型角度(相关背景知识可参考Milleretal.，2018；Ratcliff&Smith，2004)解释JILMEM的构建，如图1所示。从图1可以看出，所有刺激信息的累加起点是相同的，并且设置的时间限制为1500ms，其中实线①代表被试在刺激项目上做出了正确操作反应的信息累加过程，其累加反应时服从LMEM分布；而实线②因固定时间窗口实验中刺激呈现时间限制，或者在没有固定时间窗口实验中被试提前进行了操作反应，进而导致操作错误，此时被试的信息累加过程并不完整，理论上还需要经过虚线③的信息累加过程，才能做出正确的操作反应，其信息累加的反应时是截尾分布(Miller，1998；Milleretal.，2018)，采用截尾的LMEM分布进行描述。LMEM只关注了实线①的正确反应时的信息累加过程的分析，而忽视了实线②和虚线③的错误反应时信息累加过程的分析，两者联合构成JILMEM模型。 图1JILMEM的信息累加扩散过程 基于此，JILMEM的数学表达式可以表示为 其中，Xijk为被试i在实验条件k上对刺激项目j上的操作反应，1为正确操作，0为错误操作；f(Tijk)为Tijk的密度函数(式1)，S(Tijk)为Tijk的生存函数，通过1减去LMEM的累积分布表示，由于LMEM残差服从正态分布，LMEM的累积分布为正态累积分布；模型包含的参数与LMEM(式1)是相同的。 2.2贝叶斯参数估计 对本研究采用的模型与方法如表1所示。从表1中可以看出，LMEM分别应用于原始数据(ORT)与均值替换(MRT)反应时，同时采用R语言lmerTest包(Brown，2021)的REML方法(RestrictedMaximumLikelihoodmethod)和rstan包(Gelmanetal.，2015)的汉密尔顿蒙特卡洛(HamiltonianMonteCarlo，HMC)算法进行估计。此外，JILMEM只采用rstan包的HMC算法进行估计，并且与F1/F2检验结果进行比较。LMEM和JILMEM待估参数的先验分布设定(附录1)以及JILMEM的详细应用stan代码(附录2)，感兴趣的读者可以通过以下网址/anonymous/QnZRUmJp进行下载。 注：p指处理水平数，N指被试量，m为刺激项目量，F1-ORT/F1-MRT指分别基于原始数据ORT和均值替换MRT的F1检验，F2-ORT/F2-MRT指分别基于ORT和MRT数据的F2检验，R-ORT指LMEM基于ORT数据采用R软件lmerTest包的REML方法估计；R-MRT指LMEM基于MRT数据采用R工具包lmerTest的REML方法估计；stan-ORT指LMEM基于ORT数据采用stan的HMC方法估计；stan-MRT指LMEM基于MRT数据采用stan的HMC方法估计。下同。 3实验数据分析 3.1实验描述 本研究选用Adelman等(2014)的掩蔽形式启动的行为数据库的部分数据作为本研究的实例数据(数据来源https://files.warwick.ac.uk/jadelman2/browse#FPP)。为了获得相同被试在相同刺激类型启动的数据，本研究经过筛选最终得到36个被试在两个启动类型各30个项目刺激的数据，实验设计为单因素两水平重复测量实验设计。两种启动类型分别为启动词与目标词相同，以及启动词与目标词无任何相关。实验流程是每次试验都是以300ms的注视点(+)开始呈现，然后是500ms的(#####)掩蔽刺激呈现。之后，启动词以小写字母的五分之八大小呈现50ms，然后目标词以大写字母呈现，直到被试对于非单词或单词做出左或右反应，或经过2000ms为止。 3.2分析 F1/F2检验和LMEM分别应用原始反应时数据(originalreactiontime，ORT)和均值替换的反应时数据(mean-replacementreactiontime，MRT)；而JILMEM则联合分析正确与错误反应时数据。在HMC估计中设定2条马尔科夫链，每条链包含40000次迭代，删除前面的30000次，最后保留两条链的剩余10000次，共20000次迭代进行参数估计推断。采用潜在量尺缩减因子(potentialscalereductionfactor，PSRF)(Brooks&Gelman，1998)评估参数估计的收敛，本研究各参数的PSRF均要求小于1.1以表示估计已收敛。 在模型比较上，R语言包lmerTEST使用AIC和BIC评价模型的效果；而stan基于贝叶斯思想构建了WAIC(widelyavailableinformationcriterion)和LOO(leave-one-outcross-validation)两个模型拟合统计量(Vehtarietal.，2017)来进行评价。越小的AIC、BIC、WAIC和LOO表明模型对数据拟合越好。 3.3结果 表2呈现了不同模型的拟合结果。在不同数据上，MRT数据的LMEM模型的拟合指标AIC、BIC、WAIC和LOO都比ORT数据的LMEM模型的拟合指数要小。然而，相对于LMEM模型，联合错误反应时的JILMEM模型的拟合指数最小。结果表明，联合错误反应时的JILMEM能显著提高模型对数据的拟合效果。 表3呈现了LMEM和JILMEM的参数估计值及其置信区间。在固定效应参数β0和β1上，LMEM估计值都要小于JILMEM的参数估计值；在随机效应上，除了τ11外，LMEM随机效应参数也都要低于JILMEM。对于ORT和MRT数据估计结果，发现除了τ11和β1两个参数，基于MRT的LMEM的参数估计结果都要小于ORT的估计结果。在不同估计方法中，REML和HMC对LMEM参数估计结果接近。此外，超参数被试水平随机截距和斜率间相关ρs都包含0，且区间范围较大，该参数在已有研究也发现很难准确估计(Sorensenetal.，2016)。 4模拟研究 4.1数据生成与分析 在模拟研究中，参考实证研究估计值进行了6个条件的模拟比较(如表4所示)，此外不同效应量(d=0，0.2，0.5)与6个条件交叉，每个条件重复100次。为了使得生成的数据符合截尾数据特点，本研究设置的数据生成过程如下：在不同条件C1-C6下，首先根据I0j～N(0，ω200)、ξijk～N(0，σ2)和式(2)生成随机效应参数值，同时根据对应效应量d计算固定效应参数β1(Juddetal.，2017)以及加上固定效应β0代入式(1)得到被试在刺激项目上的正确反应时RT。然后，固定每个刺激项目呈现时间为1500ms，并基于正态分布N(1000，σ2)生成截尾时间CT。最后，当RT＜CT且RT＜1500ms时，则被试在刺激项目上的反应正确，反应时为RT；当RT＞CT且CT＜1500ms时，则被试在刺激项目上反应错误，反应时为截尾反应时CT；当RT和CT均大于1500ms，反应错误，反应时为1500ms。 4.2结果 4.2.1模型参数估计精度比较 矩阵图可以更好地反映不同模型参数估计趋势和特点，结果如图2所示。在矩阵基于线性混合效应模型的心理实验数据建模图中，每一行的图形对应不同参数的估计均值和95％置信区间，每一列的图形对应不同模拟条件；在每一个行列交叉的图形中，横坐标为不同模型和数据对应的参数估计均值和95％置信区间。从图中可以看出，相对基准条件(C1)，随着被试(C2)和刺激项目量(C3)增加，各模型的参数标准差明显下降，结果更稳定，同时各参数的均值也更接近真值；在被试随机截距(C4)和斜率(C5)以及刺激项目随机截距(C6)变异增大后，各数据模型的参数标准差明显增大，且参数估计精度略微下降。在不同数据上，通过REML(R)与HMC(stan)估计方法的LMEM在MRT数据上的参数在不同条件上，除了被试随机斜率τ11外，其他参数都要明显小于ORT数据估计结果，但是参数的标准差更小，同时REML与HMC两种参数估计方法结果接近。相对JILMEM参数估计结果，LMEM在ORT和MRT数据上的各参数明显被低估，JILMEM参数结果更接近模拟真值。 图2不同LMEM参数估计效果 4.2.2不同统计分析方法效果比较 根据模拟的效应量不同，可以比较不同分析方法的第I类错误和统计检验力。第I类错误率通常认为越接近真值0.05越好，且在0.025～0.075之间是可以接受的范围；统计检验力是估计值显著不等于0的概率，越趋近于1越好(温忠麟，汤丹丹，顾红磊，2019)。 在无效应量d=0时，不同方法的第I类错误率的结果如图3所示。从图中可以看出，在不同模拟条件上，JILMEM都在可接受范围，并都趋近0.05。对于LMEM模型，基于ORT数据的LMEM的第I类错误率也都在可接受范围内，同时比JILMEM波动更大，但是MRT数据的LMEM在不同条件上的第I类错误率明显提高，大部分超出了可接受范围。对F1/F2检验，F1在ORT数据上的所有条件的第I类错误率都明显超出可接受范围，并且在MRT数据上明显增高；F2在ORT数据上随着被试随机斜率变异(C5)增大，第I类错误率明显超出可接受范围，其他条件则都接近可接受范围，但是在MRT数据上则全部超出可接受范围。 图3不同统计分析方法的第I类错误率(％)：无效应量d=0 在低效应量d=0.2时(图4)，JILMEM的统计检验力在刺激随机截距变异增大(C6)时最低且只有0.20，但是比ORT数据的LMEM各条件的统计检验力都高，而LMEM在MRT数据上统计检验力都明显更高。对于F1/F2检验方法，在ORT数据上，F1统计量在被试随机斜率变异增大(C5)时，统计检验力最低为0.41，而F2统计量的统计检验力与JILMEM各有优劣，F1/F2方法在MRT数据上各条件的统计检验力都有所提高。 图4不同统计分析方法的统计检验力(％)：低效应量d=0.2 在中效应量d=0.5上(图5)，F2、LMEM和JILMEM在刺激项目随机截距变异增大(C6)条件上统计检验力为0.6左右，表现最差。这可能是因为刺激项目随机截距变异增大，会导致不同水平的刺激项目之间都存在较大反应时差异。例如，汉字识别中刺激项目既包含高频字，又包含低频字，但并未将词频作为自变量考虑进去。这种情况下，不同水平的刺激项目的操作反应时波动较大，使得处理效应不易被识别，进而导致较低的统计检验力。在其他条件上，统计检验力均达到0.9及以上，而F1则始终接近1，都具有较高的统计检验力。 图5不同统计分析方法的统计检验力比较(％)：中效应量d=0.5 5总结与展望 为更准确揭示实验设计特点，本研究借鉴生存分析思想，提出联合错误反应时的线性混合效应模型(JILMEM)。在实验数据中，发现联合错误反应时的JILMEM具有明显拟合优势。在模拟研究中，REML和HMC两种估计方法在LMEM参数估计上的结果接近，同时基于MRT数据比ORT数据的估计的模型参数有更大偏差，