在统计课上，我们有一个最后的作业来处理（最后lexdec的languageR完整脚本）的数据集。任务是检查几个效应是否影响受试者的阅读时间（RT）。我们从一个简单的模型开始（“RT 是否受试验次数的影响？”）并扩展它（“NativeLanguage 有什么影响”等），直到我们最终得到这个：
model4 = lmer(RT ~ Trial + NativeLanguage + (1|Word) + (1 + Trial|Subject), lexdec)
这是正确的（因为这不是实际分析，而是一个玩具示例，我们现在可以忽略有关模型未收敛的警告......）。然而，我和我的团队合作伙伴想知道更多：如何阅读和解释输出以检查我们尝试使用模型预测的效果是否显着？
我们学到的一种方法是 t 检验。在没有混合效应的最简单模型中，我们免费获得结果：
summary(lm(RT ~ Trial, lexdec))
产生可解释的
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.4172041 0.0144549 443.948 <2e-16 ***
Trial -0.0003060 0.0001256 -2.435 0.015 *
甜的！如果我们将 alpha 值固定为标准 0.05（或 0.025，因为它是一个双边测试），那么 Trial 的影响实际上是显着的。
但lmer不提供这个，而是我们只得到报告的 t 值。所以这是我们检查这些值是否重要的脚本：
ttest <- function (model, effect, df, alpha=0.05) {
# Performs a two-sided t-test for a given alpha value on the given model.
t <- abs(coef(summary(model))[effect, 't value'])
q <- qt(1 - alpha / 2, df=df)
paste('|t| =', round(t, 3), '>', round(q, 3), '?', t > q)
}
现在棘手的部分：如何提出自由度（df）？我一直将 df 理解为我可以独立排列的变量的数量（有时会限制其他变量，例如机器人关节）。
很多时候，我会被以下定义抛出：（参见此处的多重线性回归自由度），当然还有学生和韦尔奇 t 检验的 df 公式。所有这些都会导致不同数量的 df，但是对于无限多的 df，它们无论如何都无关紧要（维基百科的 t 值表很好地表明，即使对于 100+，它对于“我们的”小型行为实验来说也几乎不重要）。
话虽如此，这样的库的作者已经解决了这个问题似乎是一个合理的假设，是的，在summary（lm！）中也有许多 dfs: F-statistic: 5.931 on 1 and 1657 DF, p-value: 0.01498。等等，1657？这似乎很多。但它是样本数（）减去 2，因此：和，它必须是唯一的变量。Trial
所以我假设图书馆的样本数量是最好的猜测。但是：如果我在测试我的科目之前干净地编写我的模型，我将不知道我最终会得到多少 df（我会有 30 个科目吗？还是只有 29 个，因为一个人生病了？），因为这个数字取决于。所以我的模型中一定有一个变量，当然我可以找到同意我的人（线性混合模型中的自由度计算）：
在 [...] 模型中，您有一个截距参数、每个随机效应项的方差参数和给定随机效应的观测值的方差参数
如果我数一下，我会得出 7 个：截距、试验方差、NativeLanguage-Variance、单词截距-方差、主题-截距-方差、主题-试验-方差、残差。
这看起来更清晰，但不幸的是，它改变了我们很多重要的结果，即使对于简单的模型也是如此。虽然内置lm告诉我们：Pr(>|t|) = 0.015，但我们现在找到：|t| = 2.435 > 4.303 ? FALSE，因为 df 现在太小了。
此外，对于复杂模型，我们得到多个 t 值。我们应该考虑哪些？如果一切都很重要，那么重要的事情吗？一？我听说了很多关于多项测试的更正 - 如果我测试多个参数，是否已经需要它？
有人可以阐明我陷入的这种意义和黑暗吗？我的方法甚至正确吗？或者没有简单的答案，我总是必须考虑许多不同的东西作为 df - 并且可以通过使用足够的样本和参与者来假设无限多（我有时听说拥有太多也不是很好，没有解释，但我认为这是由于过度拟合？）？
一个类似的问题是如何确定序数回归中 t 检验的自由度？，没有答案，我很享受How to get the degree of Freedom when doing multiple linear regression? ，不幸的是，它比较了嵌套模型——我目前对此不感兴趣。
这是一个完整的工作示例代码，看看我在说什么，然后是它的输出：
library('languageR')
library('lme4')
ttest <- function (model, effect, df, alpha=0.05) {
# Performs a two-sided t-test for a given alpha value on the given model.
t = abs(coef(summary(model))[effect, 't value'])
q = qt(1 - alpha / 2, df=df)
paste('|t| =', round(t, 3), '>', round(q, 3), '?', t > q)
}
model1 = lm(RT ~ Trial, lexdec)
summary(model1)
print(paste("Is reaction time affected by the number of trials only?", ttest(model1, 'Trial', 2)))
# snip models 2-3
model4 = lmer(RT ~ Trial + (1 + Trial|Subject) + (1|Word) + NativeLanguage, lexdec)
summary(model4)
print(paste("Is reaction time still affected by trials if considering native language and by-subject variance?", ttest(model4, 'Trial', 7)))
print(paste("Is reaction time still affected by native language if considering by-subject variance?", ttest(model4, 'NativeLanguageOther', 7)))
输出：
> library('languageR')
> library('lme4')
>
> ttest <- function (model, effect, df, alpha=0.05) {
+ # Performs a two-sided t-test for a given alpha value on the given model.
+ t = abs(coef(summary(model))[effect, 't value'])
+ q = qt(1 - alpha / 2, df=df)
+ paste('|t| =', round(t, 3), '>', round(q, 3), '?', t > q)
+ }
>
> model1 = lm(RT ~ Trial, lexdec)
> summary(model1)
Call:
lm(formula = RT ~ Trial, data = lexdec)
Residuals:
Min 1Q Median 3Q Max
-0.53961 -0.16943 -0.04366 0.11594 1.18357
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.4172041 0.0144549 443.948 <2e-16 ***
Trial -0.0003060 0.0001256 -2.435 0.015 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2412 on 1657 degrees of freedom
Multiple R-squared: 0.003567, Adjusted R-squared: 0.002965
F-statistic: 5.931 on 1 and 1657 DF, p-value: 0.01498
> print(paste("Is reaction time affected by the number of trials only?", ttest(model1, 'Trial', 2)))
[1] "Is reaction time affected by the number of trials only? |t| = 2.435 > 4.303 ? FALSE"
>
> # snip models 2-3
>
> model4 = lmer(RT ~ Trial + (1 + Trial|Subject) + (1|Word) +NativeLanguage, lexdec)
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00759802 (tol = 0.002, component 1)
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?
> summary(model4)
Linear mixed model fit by REML ['lmerMod']
Formula: RT ~ Trial + (1 + Trial | Subject) + (1 | Word) + NativeLanguage
Data: lexdec
REML criterion at convergence: -922.6
Scaled residuals:
Min 1Q Median 3Q Max
-2.3256 -0.6159 -0.0954 0.4647 6.0361
Random effects:
Groups Name Variance Std.Dev. Corr
Word (Intercept) 5.942e-03 0.0770873
Subject (Intercept) 3.179e-02 0.1783087
Trial 4.937e-07 0.0007026 -0.72
Residual 2.869e-02 0.1693814
Number of obs: 1659, groups: Word, 79; Subject, 21
Fixed effects:
Estimate Std. Error t value
(Intercept) 6.3403360 0.0477379 132.82
Trial -0.0002435 0.0001778 -1.37
NativeLanguageOther 0.1640987 0.0563356 2.91
Correlation of Fixed Effects:
(Intr) Trial
Trial -0.609
NtvLnggOthr -0.506 0.001
convergence code: 0
Model failed to converge with max|grad| = 0.00759802 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?
> print(paste("Is reaction time still affected by trials if considering native language and by-subject variance?", ttest(model4, 'Trial', 7)))
[1] "Is reaction time still affected by trials if considering native language and by-subject variance? |t| = 1.369 > 2.365 ? FALSE"
> print(paste("Is reaction time still affected by native language if considering by-subject variance?", ttest(model4, 'NativeLanguageOther', 7)))
[1] "Is reaction time still affected by native language if considering by-subject variance? |t| = 2.913 > 2.365 ? TRUE"