混合模型中的计划对比

Question

我的实验由两个因素组成，一个是主体间的因素，另一个是学科内的因素。Age是主体间因素，有两个层次：(低和高)。Time是主体内的因素，包含三个层次: 1，4，5。dv是我的因变量，id是每个参与者的标识符。我附上了6名第一批参与者的数据。

使用R，我进行了一个方差分析，对这两个因素都产生了显著的结果。我有两个计划中的对比：

1. 在第一次，我想比较这两个年龄组(低和高)。也就是说，主体之间的比较。
2. 在年龄低的时候，我想比较一下时间1和时间5。也就是说，在主题比较中。

当然，我可以执行t测试，但这似乎不合适，因为我可以将我的标准错误估计建立在更多的单元格上。我的问题是，我如何才能做到上述的对比，什么是适当的自由程度？

structure(list(id = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 
5L, 5L, 5L, 5L, 5L, 5L), .Label = c("1", "2", "3", "5", "6", 
"7", "8", "11", "12", "13", "15", "17", "18", "19", "20", "21", 
"22", "23", "24", "25", "27", "28", "29", "31", "32", "34", "35", 
"36", "37", "38", "39", "40", "42", "43", "44", "45", "46", "47", 
"48", "49", "52", "53", "54", "55", "56", "58", "59", "60", "62", 
"63", "64", "66", "67", "68", "69", "70", "71", "72", "73", "74", 
"75", "77", "79", "80", "81", "83", "84", "85", "86", "87", "88", 
"89", "90", "91", "92", "93", "94", "96", "97", "98", "99", "100", 
"101", "102", "103", "104", "105", "106", "107", "108", "109", 
"110", "111", "112", "113", "114", "115", "116", "117", "118", 
"119", "120", "121", "122", "123", "124", "125", "126", "127", 
"128", "129", "130", "132", "133", "134", "135", "136", "137", 
"138", "139", "140", "142", "143", "144", "145", "146", "147", 
"148", "149", "150", "151", "152", "153", "154", "156", "157", 
"158", "159", "160", "161", "162", "163", "165", "166", "167", 
"168", "169", "171", "172", "174", "175", "176", "177", "178", 
"179", "180", "181", "182", "183", "184", "185", "186", "187", 
"188", "189", "190", "191", "192", "193", "194", "195", "196", 
"200", "201", "202", "203", "204", "205", "206", "208", "209", 
"210", "212", "213", "214", "215", "216", "217", "218", "219", 
"220", "222", "223", "224", "226", "228", "230", "231", "232", 
"233", "234", "236", "237", "238", "239", "240", "241", "242", 
"243", "244", "246", "247", "248", "249", "250", "251", "252", 
"253", "254", "255", "256", "257", "258", "260", "261", "262", 
"263", "266", "267", "269", "270", "271", "272", "273", "274", 
"275", "276", "277", "278", "279", "280", "281", "282", "283", 
"284", "285", "286", "287", "288", "289", "290", "291", "292", 
"293", "294", "295", "296", "298", "299", "300"), class = "factor"), 
    age = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 2L, 2L, 2L), .Label = c("high", "low"), class = "factor"), 
    time = structure(c(3L, 3L, 2L, 2L, 1L, 1L, 3L, 3L, 2L, 2L, 
    1L, 1L, 3L, 3L, 2L, 2L, 1L, 1L, 3L, 3L, 2L, 2L, 1L, 1L, 3L, 
    3L, 2L, 2L, 1L, 1L), .Label = c("1", "4", "5"), class = "factor"), 
    dv = c(104, 102, 104, 103, 104, 104, 102, 102, 102, 102, 
    106, 106, 106, 106, 107, 107, 106, 106, 106, 107, 105, 104, 
    106, 107, 104, 101, 104, 101, 104, 106)), row.names = c(NA, 
-30L), class = c("tbl_df", "tbl", "data.frame"))

Answer 1

加载适合多级模型的lme4包(重复度量方差分析的替代方案)和emmeans，后者可以进行对比。

library(lme4)
library(emmeans)

这符合一个模型，其中dv是由交互(R自动填充-在主要影响)，加上一个随机拦截和time-both嵌套在id中的随机效应预测的。dat是我在你的帖子中从dput中保存下来的东西。

mod <- lmer(dv ~ age * time + (1 + time | id), dat)

反差令人困惑，我总是偏执狂--我会弄错的。所以我们可以用emmeans找到它们。我们可以拟合一个emmeans对象，它可以在time和age的每个组合中获取值。

emm_mod <- emmeans(mod, ~ time + age)

我们想要的对比是pairs()结果的第三个和第十四个(自己运行它，看看它看起来是什么样子)。通过将coef()放在pairs()对象周围，可以得到所需的特定对比。你只想要两列-第三和第十四栏：

(contr_mat <- coef(pairs(emm_mod))[, c("c.3", "c.14")])

返回：

       c.3 c.14
1,high   1    0
4,high   0    0
5,high   0    0
1,low   -1    1
4,low    0    0
5,low    0   -1

您可以通过在contr中指定它们来关注这两个对比。您也可以使用adjust选择的p值调整到任何东西--我将从"holm"开始。

emmeans(mod, ~ time + age, contr = contr_mat, adjust = "holm")

contrasts位将提供您感兴趣的p值：

$emmeans
 time age    emmean        SE df lower.CL upper.CL
 1    high 106.2500 0.6518719  3 104.1755 108.3245
 4    high 105.7500 0.8544406  3 103.0308 108.4692
 5    high 106.2500 0.4759431  3 104.7353 107.7647
 1    low  105.0000 0.5322511  3 103.3061 106.6939
 4    low  102.6667 0.6976478  3 100.4464 104.8869
 5    low  102.5000 0.3886059  3 101.2633 103.7367

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
 contrast estimate        SE df t.ratio p.value
 c.3          1.25 0.8415630  3   1.485  0.2341
 c.14         2.50 0.7104068  3   3.519  0.0779

P value adjustment: holm method for 2 tests 

您也可以尝试adjust = "none"

> emmeans(mod, ~ time + age, contr = contr_mat, adjust = "none")$contrasts
 contrast estimate        SE df t.ratio p.value
 c.3          1.25 0.8415630  3   1.485  0.2341
 c.14         2.50 0.7104068  3   3.519  0.0389