r - glmer 中的 glht：需要重新调平吗？
问题描述
我使用 glmer 来分析重复测量设计中的反应时间，其中包含 2 个固定因素（每个 3 个水平）和 1 个随机因素。由于我有一个重要的互动，我想做事后得到对比。我首先使用 relevel 来根据需要多次更改模型的参考，以获得我需要的对比度。我已经看到 glht 也可以使用（并避免重新调整），但我得到一些在查看数据时没有任何意义的 p 值。我尝试在 glht 之前重新调整我的一个因素，并且 p 值更有意义。我的对比度矩阵好吗？这是否意味着我还需要在 glht 中使用 relevel？任何人都可以帮助我吗？
我有一个实验，有 2 个重复的固定因子 Verb (vt,vnt,vnta) 和 Delay (170,350,500)，而 Subjects 作为随机因子。在测量反应时间时，我按照 Lo & Andrews (2015) 的建议使用了 glmer（lme4 包）。模型的方差分析给出了 2 个固定因素之间的显着交互作用，我会进行事后测试。我首先使用 relevel 来更改模型的参考并测试所有需要的对比。我发现某些条件之间存在显着差异，但是我不确定是否对数据进行了任何更正。这是我的第一个问题。
我已经看到 glht（multcomp 包）也可以用于事后处理（并避免重新调整）。因此，我尝试了这个，但我不确定我是否做得很好，因为在查看数据时我得到了看起来非常奇怪的显着差异。当我在 glht 之前使用 relevel 时，我得到了一些更连贯的结果（以及与第一次 relevel 分析中类似的估计）。这是否意味着 glht 有时需要使用 relevel ？谁能帮我？
> m2b <- glmer(RT ~ Delai * Verbe + (1|Sujet), data = Data_Tacti2, family=Gamma(link="identity"), glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
> Anova(m2b)
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: RT
Chisq Df Pr(>Chisq)
Delai 677.379 2 <2e-16 ***
Verbe 3.562 2 0.1685
Delai:Verbe 10.161 4 0.0378 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
型号参考 VNT 170
> summary(m2b)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Gamma ( identity )
Formula: RT ~ Delai * Verbe + (1 | Sujet)
Data: Data_Tacti2
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
AIC BIC logLik deviance df.resid
51151.3 51220.9 -25564.6 51129.3 4122
Scaled residuals:
Min 1Q Median 3Q Max
-3.3223 -0.6339 -0.1768 0.4510 6.8590
Random effects:
Groups Name Variance Std.Dev.
Sujet (Intercept) 836.62309 28.9244
Residual 0.08911 0.2985
Number of obs: 4133, groups: Sujet, 18
Fixed effects:
Estimate Std. Error t value Pr(>|z|)
(Intercept) 494.110 9.414 52.485 < 2e-16 ***
Delai350 -76.068 4.328 -17.576 < 2e-16 ***
Delai500 -93.102 4.313 -21.585 < 2e-16 ***
Verbevnta -7.559 4.901 -1.542 0.12298
Verbevt -16.837 5.224 -3.223 0.00127 **
Delai350:Verbevnta 7.926 6.657 1.191 0.23379
Delai500:Verbevnta 8.381 6.271 1.337 0.18136
Delai350:Verbevt 19.670 7.014 2.804 0.00504 **
Delai500:Verbevt 10.198 5.624 1.813 0.06980 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Del350 Del500 Vrbvnt Verbvt Dl350:Vrbvn Dl500:Vrbvn Dl350:Vrbvt
Delai350 -0.262
Delai500 -0.212 0.381
Verbevnta -0.232 0.206 0.139
Verbevt -0.297 0.196 0.168 0.184
Dl350:Vrbvn 0.209 -0.360 -0.019 -0.555 0.024
Dl500:Vrbvn 0.280 -0.182 -0.439 -0.515 -0.051 0.351
Dl350:Vrbvt 0.297 -0.408 -0.066 -0.011 -0.576 0.101 0.038
Dl500:Vrbvt 0.195 -0.083 -0.350 0.052 -0.571 -0.106 0.131 0.406
重新调平至 VT 170
> Data_Tacti2$Delai <- relevel (Data_Tacti2$Delai, ref = "170")
> Data_Tacti2$Verbe <- relevel (Data_Tacti2$Verbe, ref = "vt")
> m2b_vt_170 <- glmer(RT ~ Delai * Verbe + (1|Sujet), data = Data_Tacti2, family=Gamma(link="identity"), glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
> summary(m2b_vt_170)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Gamma ( identity )
Formula: RT ~ Delai * Verbe + (1 | Sujet)
Data: Data_Tacti2
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
AIC BIC logLik deviance df.resid
51151.3 51220.9 -25564.6 51129.3 4122
Scaled residuals:
Min 1Q Median 3Q Max
-3.3223 -0.6339 -0.1768 0.4510 6.8590
Random effects:
Groups Name Variance Std.Dev.
Sujet (Intercept) 836.63773 28.9247
Residual 0.08911 0.2985
Number of obs: 4133, groups: Sujet, 18
Fixed effects:
Estimate Std. Error t value Pr(>|z|)
(Intercept) 477.272 7.337 65.050 < 2e-16 ***
Delai350 -56.397 4.444 -12.692 < 2e-16 ***
Delai500 -82.904 4.670 -17.753 < 2e-16 ***
Verbevnt 16.838 4.189 4.019 5.83e-05 ***
Verbevnta 9.279 5.293 1.753 0.07962 .
Delai350:Verbevnt -19.671 5.984 -3.287 0.00101 **
Delai500:Verbevnt -10.198 7.293 -1.398 0.16200
Delai350:Verbevnta -11.745 6.581 -1.785 0.07433 .
Delai500:Verbevnta -1.817 5.805 -0.313 0.75428
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Del350 Del500 Verbevnt Verbevnta Delai350:Verbevnt Delai500:Verbevnt Delai350:Verbevnta
Delai350 -0.114
Delai500 -0.015 0.276
Verbevnt 0.012 -0.010 0.108
Verbevnta -0.157 0.220 0.070 0.139
Delai350:Verbevnt 0.006 -0.226 0.095 -0.420 0.091
Delai500:Verbevnt 0.003 0.107 -0.456 -0.438 0.152 0.183
Delai350:Verbevnta 0.198 -0.428 0.006 0.101 -0.570 0.007 -0.177
Delai500:Verbevnta 0.062 -0.038 -0.332 0.043 -0.511 -0.173 0.100 0.342
重新调平至 VNT 350
> Data_Tacti2$Delai <- relevel (Data_Tacti2$Delai, ref = "350")
> Data_Tacti2$Verbe <- relevel (Data_Tacti2$Verbe, ref = "vnt")
> m2b_vnt_350 <- glmer(RT ~ Delai * Verbe + (1|Sujet), data = Data_Tacti2, family=Gamma(link="identity"), glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
> summary(m2b_vnt_350)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Gamma ( identity )
Formula: RT ~ Delai * Verbe + (1 | Sujet)
Data: Data_Tacti2
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
AIC BIC logLik deviance df.resid
51151.3 51220.9 -25564.6 51129.3 4122
Scaled residuals:
Min 1Q Median 3Q Max
-3.3223 -0.6339 -0.1768 0.4510 6.8590
Random effects:
Groups Name Variance Std.Dev.
Sujet (Intercept) 836.62212 28.9244
Residual 0.08911 0.2985
Number of obs: 4133, groups: Sujet, 18
Fixed effects:
Estimate Std. Error t value Pr(>|z|)
(Intercept) 420.875 7.076 59.478 < 2e-16 ***
Delai170 56.398 5.221 10.802 < 2e-16 ***
Delai500 -26.507 4.087 -6.486 8.8e-11 ***
Verbevnt -2.833 4.184 -0.677 0.498329
Verbevnta -2.466 4.030 -0.612 0.540484
Delai170:Verbevnt 19.671 5.219 3.769 0.000164 ***
Delai500:Verbevnt 9.473 5.672 1.670 0.094900 .
Delai170:Verbevnta 11.745 5.454 2.153 0.031280 *
Delai500:Verbevnta 9.928 5.097 1.948 0.051439 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Del170 Del500 Verbevnt Verbevnta Delai170:Verbevnt Delai500:Verbevnt Delai170:Verbevnta
Delai170 0.032
Delai500 0.107 0.225
Verbevnt -0.022 -0.140 0.117
Verbevnta -0.010 -0.035 0.073 0.235
Delai170:Verbevnt 0.041 -0.310 -0.045 -0.229 0.066
Delai500:Verbevnt -0.096 0.200 -0.356 -0.468 0.029 0.137
Delai170:Verbevnta -0.087 -0.438 0.037 0.181 -0.115 0.115 -0.189
Delai500:Verbevnta -0.105 0.094 -0.336 -0.045 -0.383 -0.039 0.184 -0.019
我根据需要多次使用 relevel 来获得因子动词的 3 个级别和因子延迟的 3 个级别之间的所有对比（此处未报告）
现在使用 glht
型号 m2b，参考 VT 170
> contrast.matrix <- rbind(
+ `VNT170 vs VT170` = c(0, 0, 0, 0, 1, 0, 0, 0, 0),
+ `VNTA170 vs VT170` = c(0, 0, 0, -1, 1, 0, 0, 0, 0),
+ `VNT170 vs VNTA170` = c(0, 0, 0, 1, 0, 0, 0, 0, 0),
+ `VNT350 vs VT350` = c(0, 1, 0, 0, 0, 0, 0, -1, 0),
+ `VNTA350 vs VT350` = c(0, 0, 0, 0, 0, 1, 0, -1, 0),
+ `VNT350 vs VNTA350` = c(0, 1, 0, 0, 0, -1, 0, 0, 0),
+ `VNT500 vs VT500` = c(0, 0, 1, 0, 0, 0, 0, 0, -1),
+ `VNTA500 vs VT500` = c(0, 0, 0, 0, 0, 0, 1, 0, -1),
+ `VNT500 vs VNTA500` = c(0, 0, 1, 0, 0, 0, -1, 0, 0)
+ )
> comps_refVNT170 <- glht(m2b, contrast.matrix)
> summary(comps_refVNT170)
Simultaneous Tests for General Linear Hypotheses
Fit: glmer(formula = RT ~ Delai * Verbe + (1 | Sujet), data = Data_Tacti2,
family = Gamma(link = "identity"), control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 1e+05)))
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
VNT170 vs VT170 == 0 -16.837 5.224 -3.223 0.010 *
VNTA170 vs VT170 == 0 -9.278 6.472 -1.434 0.621
VNT170 vs VNTA170 == 0 -7.559 4.901 -1.542 0.544
VNT350 vs VT350 == 0 -95.738 9.629 -9.943 <0.001 ***
VNTA350 vs VT350 == 0 -11.744 9.171 -1.281 0.726
VNT350 vs VNTA350 == 0 -83.994 9.152 -9.178 <0.001 ***
VNT500 vs VT500 == 0 -103.300 8.198 -12.601 <0.001 ***
VNTA500 vs VT500 == 0 -1.816 7.856 -0.231 1.000
VNT500 vs VNTA500 == 0 -101.484 9.038 -11.228 <0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
我得到对比 VNT350-VT350 的 p 值 <.001，这非常奇怪（查看数据时 VNT500-VT500 和 VNT500-VNTA500 相同），并且这种差异在第一次重新调平分析中并不显着
矩阵好吗？
如果我重新调整对 VNT 350 的 tp 参考（尽管我不应该这样做），则估计与第一次重新调整分析中的估计相似，并且 VNT350-VT350 的 p 值更有意义
这是否意味着我也需要为 glht 重新调平？
> Data_Tacti2$Delai <- relevel (Data_Tacti2$Delai, ref = "350")
> Data_Tacti2$Verbe <- relevel (Data_Tacti2$Verbe, ref = "vnt")
> m2b_vnt_350 <- glmer(RT ~ Delai * Verbe + (1|Sujet), data = Data_Tacti2, family=Gamma(link="identity"), glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
> summary(m2b_vnt_350)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Gamma ( identity )
Formula: RT ~ Delai * Verbe + (1 | Sujet)
Data: Data_Tacti2
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
AIC BIC logLik deviance df.resid
51151.3 51220.9 -25564.6 51129.3 4122
Scaled residuals:
Min 1Q Median 3Q Max
-3.3223 -0.6339 -0.1768 0.4510 6.8590
Random effects:
Groups Name Variance Std.Dev.
Sujet (Intercept) 836.62212 28.9244
Residual 0.08911 0.2985
Number of obs: 4133, groups: Sujet, 18
Fixed effects:
Estimate Std. Error t value Pr(>|z|)
(Intercept) 420.875 7.076 59.478 < 2e-16 ***
Delai170 56.398 5.221 10.802 < 2e-16 ***
Delai500 -26.507 4.087 -6.486 8.8e-11 ***
Verbevnt -2.833 4.184 -0.677 0.498329
Verbevnta -2.466 4.030 -0.612 0.540484
Delai170:Verbevnt 19.671 5.219 3.769 0.000164 ***
Delai500:Verbevnt 9.473 5.672 1.670 0.094900 .
Delai170:Verbevnta 11.745 5.454 2.153 0.031280 *
Delai500:Verbevnta 9.928 5.097 1.948 0.051439 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Del170 Del500 Verbevnt Verbevnta Delai170:Verbevnt Delai500:Verbevnt Delai170:Verbevnta
Delai170 0.032
Delai500 0.107 0.225
Verbevnt -0.022 -0.140 0.117
Verbevnta -0.010 -0.035 0.073 0.235
Delai170:Verbevnt 0.041 -0.310 -0.045 -0.229 0.066
Delai500:Verbevnt -0.096 0.200 -0.356 -0.468 0.029 0.137
Delai170:Verbevnta -0.087 -0.438 0.037 0.181 -0.115 0.115 -0.189
Delai500:Verbevnta -0.105 0.094 -0.336 -0.045 -0.383 -0.039 0.184 -0.019
> contrast.matrix2 <- rbind(
+ `VNT170 vs VT170` = c(0, 1, 0, 0, 0, 0, 0, -1, 0),
+ `VNTA170 vs VT170` = c(0, 0, 0, 0, 0, 1, 0, -1, 0),
+ `VNT170 vs VNTA170` = c(0, 1, 0, 0, 0, -1, 0, 0, 0),
+ `VNT350 vs VT350` = c(0, 0, 0, 0, 1, 0, 0, 0, 0),
+ `VNTA350 vs VT350` = c(0, 0, 0, 1, -1, 0, 0, 0, 0),
+ `VNT350 vs VNTA350` = c(0, 0, 0, 1, 0, 0, 0, 0, 0),
+ `VNT500 vs VT500` = c(0, 0, 1, 0, 0, 0, 0, 0, -1),
+ `VNTA500 vs VT500` = c(0, 0, 0, 0, 0, 0, 1, 0, -1),
+ `VNT500 vs VNTA500` = c(0, 0, 1, 0, 0, 0, -1, 0, 0)
+ )
> comps_refVNT350 <- glht(m2b_vnt_350, contrast.matrix2)
> summary(comps_refVNT350)
Simultaneous Tests for General Linear Hypotheses
Fit: glmer(formula = RT ~ Delai * Verbe + (1 | Sujet), data = Data_Tacti2,
family = Gamma(link = "identity"), control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 1e+05)))
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
VNT170 vs VT170 == 0 44.6531 9.0536 4.932 < 1e-04 ***
VNTA170 vs VT170 == 0 7.9262 7.0999 1.116 0.858353
VNT170 vs VNTA170 == 0 36.7269 8.4492 4.347 0.000132 ***
VNT350 vs VT350 == 0 -2.4665 4.0296 -0.612 0.991632
VNTA350 vs VT350 == 0 -0.3665 5.0804 -0.072 1.000000
VNT350 vs VNTA350 == 0 -2.8330 4.1839 -0.677 0.985839
VNT500 vs VT500 == 0 -36.4355 7.5282 -4.840 < 1e-04 ***
VNTA500 vs VT500 == 0 -0.4557 6.8929 -0.066 1.000000
VNT500 vs VNTA500 == 0 -35.9798 8.0848 -4.450 < 1e-04 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
解决方案
推荐阅读
- mongodb - 如何重命名 MongoDB 数据库？（当前状态）
- powershell - 从日志文件中获取最后 30 分钟
- javascript - java脚本websocket无法连接到python server
- json - Json Response Spring boot 中的继承
- python - 如何使用请求模块从 div 获取文本？
- javascript - NodeJS Firebase 等待请求
- azure-devops - SAPUI5 xhr 请求 Azure DevOps 授权失败
- php - PHP Scraper 循环直到超过时间限制（使用简单的 HTML Dom 解析器）
- r - Dredge MuMIn：在 GLMM 上使用挖泥机时，空模型是否包含随机效应？
- macos - OSX上的“dyld：警告，LC_RPATH @executable_path”是什么意思，我该如何解决？