评估嵌套效应的成对比较
Evaluating pairwise comparisons of nested effects
我有一个模型，我正在尝试研究嵌套效应的某些成对比较。我不确定我是否正确编写了模型，而且我不了解如何实际评估嵌套项。
我的数据框有一个名为 'quality' 的响应变量和三个名为 "site" "month" 和 "day" 的预测变量。在我的实验设置中，我测量了每个人的质量。有两个站点。我在 4 个月内对每个网站进行了抽样。每个月有连续 4 天的抽样。我想知道某一天的个体是否与同一个月内其他日子的个体具有显着不同的质量。我对比较一个月的天数和另一个月的天数不感兴趣。
我的数据框如下
test<- structure(list(Site = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("H", "W"), class = "factor"),
Day = structure(c(19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 15L, 15L, 15L,
15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L,
15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 26L, 26L, 26L, 26L,
26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L,
26L, 26L, 8L, 8L, 8L, 8L, 8L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L,
17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L,
17L, 17L, 17L, 17L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L,
14L, 14L, 14L, 14L, 14L, 25L, 25L, 25L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 22L,
22L, 7L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,
16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,
16L, 16L, 16L, 16L, 16L, 16L, 13L, 13L, 13L, 13L, 13L, 13L,
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 24L, 24L, 24L,
24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L,
21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L,
21L, 21L, 21L, 21L, 21L, 21L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L,
23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L,
23L, 23L, 23L, 23L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L,
20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L,
20L, 20L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L), .Label = c("H1", "H10", "H11", "H12", "H13", "H14",
"H15", "H16", "H2", "H3", "H4", "H6", "H7", "H8", "H9", "W10",
"W11", "W12", "W13", "W2", "W3", "W4", "W6", "W7", "W8",
"W9"), class = "factor"), Month = structure(c(3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 2L, 2L, 2L, 2L, 2L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("August",
"November", "October", "September"), class = "factor"), Quality = c(42.535,
46.651, 45.466, 43.483, 44.896, 46.581, 47.494, 47.529, 46.562,
45.111, 45.982, 48.367, 47.39, 45.388, 46.313, 44.732, 48.641,
46.614, 45.234, 45.96, 44.795, 44.333, 46.559, 46.826, 44.166,
45.19, 46.661, 45.481, 46.828, 43.487, 49.505, 48.558, 45.218,
44.802, 43.975, 47.23, 44.85, 46.213, 44.726, 43.036, 47.211,
45.536, 44.62, 44.297, 36.115, 39.314, 42.349, 44.919, 46.296,
46.317, 45.858, 45.036, 45.861, 48.85, 45.337, 45.03, 47.4,
48.78, 49.829, 45.12, 45.599, 43.235, 44.735, 44.889, 45.666,
46.475, 44.888, 46.215, 42.242, 46.341, 45.992, 43.549, 46.612,
44.232, 42.706, 42.064, 43.837, 43.351, 41.064, 44.364, 42.597,
45.561, 44.51, 45.184, 44.896, 45.772, 47.43, 44.08, 44.697,
45.141, 43.776, 47.175, 46.115, 43.39, 47.426, 47.636, 43.672,
45.987, 45.338, 46.644, 42.192, 47.011, 45.856, 44.764, 36.285,
33.741, 34.324, 35.101, 46.844, 42.52, 48.649, 44.364, 44.688,
45.822, 44.945, 44.311, 44.684, 42.787, 45.516, 46.16, 46.289,
45.661, 45.772, 43.845, 48.717, 46.567, 44.719, 46.585, 45.33,
45.995, 48.053, 44.734, 51.233, 44.597, 45.742, 46.567, 46.478,
44.382, 47.316, 46.205, 45.111, 47.575, 46.014, 44.533, 45.347,
45.983, 47.053, 44.855, 48.021, 45.155, 49.248, 45.634, 48.815,
45.413, 43.091, 47.854, 45.19, 47.495, 47.323, 48.076, 44.183,
43.182, 46.267, 41.58, 44.237, 45.607, 48.517, 44.639, 44.773,
42.787, 43.965, 46.629, 46.256, 47.688, 44.126, 44.712, 47.097,
44.561, 47.306, 45.323, 46.328, 45.832, 46.075, 46.778, 47.445,
45.582, 47.691, 45.193, 48.453, 46.301, 44.847, 43.675, 46.066,
47.896, 45.2, 44.959, 47.401, 46.267, 45.743, 47.411, 46.926,
46.24, 46.212, 44.988, 36.552, 38.027, 47.355, 40.147, 38.094,
39.043, 37.589, 46.491, 46.413, 43.92, 45.228, 46.319, 44.764,
47.376, 43.924, 45.203, 45.418, 45.684, 46.34, 43.655, 44.365,
46.927, 48.269, 45.473, 46.451, 42.752, 48.346, 47.832, 46.534,
46.47, 43.282, 47.749, 44.856, 46.551, 45.925, 45.669, 47.263,
44.367, 47.017, 42.922, 44.904, 48.85, 45.535, 48.512, 46.154,
47.306, 46.571, 46.619, 46.092, 43.808, 47.7, 48.482, 44.407,
45.442, 44.771, 46.373, 47.777, 43.012, 46.154, 45.203, 46.443,
43.461, 45.714, 40.776, 48.949, 45.72, 48.269, 45.782, 43.945,
45.382, 43.729, 44.187, 45.267, 46.012, 42.234, 43.431, 41.973,
45.597, 45.993, 46.303, 44.493, 44.981, 46.487, 45.01, 47.009,
46.904, 48.277, 48.585, 48.625, 47.511, 44.011, 42.21, 47.124,
44.244, 47.76, 47.299, 45.278, 45.564, 44.621, 46.75, 45.396,
44.947, 46.185, 45.399, 46.095, 49.545, 47.211, 43.613, 48.494,
44.102, 45.888, 45.473, 47.222, 46.681, 45.863, 47.834, 48.386,
46.979, 46.318, 46.061, 46.347, 47.976, 47.079, 48.254, 47.643,
46.244, 46.717, 44.574, 45.177, 44.879, 46.485, 47.416, 50.235,
45.626, 48.117, 44.529, 44.281, 47.087, 47.356, 43.234, 45.841,
43.487, 42.997, 35.322, 45.554, 44.973, 43.396, 43.023, 44.65,
47.088, 41.934, 45.704, 44.559, 37.969, 42.687, 42.995, 45.287,
45.21, 43.335, 46.892, 45.534, 44.19, 43.606, 44.173, 49.334,
44.888, 47.477, 47.054, 41.041, 46.629, 45.049, 44.478, 40.278,
43.044, 43.575, 46.194, 42.688, 41.361, 46.828, 45.534, 47.395,
45.431, 45.433, 45.331, 43.947, 47.371, 48.308, 45.726, 41.833,
45.782, 44.756, 45.406, 45.661, 43.447, 46.932, 45.495, 44.349,
40.493, 43.603, 48.151, 44.037, 44.379, 45.934, 44.854, 42.321,
46.198, 44.622, 46.077, 45.306, 48.951, 47.972, 42.581, 43.608,
45.988, 44.955, 45.097, 46.768, 44.722, 45.971, 46.612, 48.956,
47.669, 47.757, 47.189, 44.184, 48.464, 49.546, 48.021, 45.448,
45.573, 46.778, 45.769, 45.419, 45.277, 47.489, 46.762, 46.238,
47.509, 47.249, 46.243, 46.124, 46.801, 47.385, 43.614, 44.661,
45.96, 48.791, 47.872, 42.402, 45.651, 45.927, 43.781, 49.923,
47.153, 46.87, 43.767, 47.3, 46.897, 44.932, 45.135, 50.124,
45.366, 45.063, 45.958, 46.731, 43.863, 45.095, 47.755, 45.446,
45.145, 45.998, 46.377, 44.369, 46.485, 48.852, 45.365, 45.934,
44.856, 48.195, 45.424, 49.05, 46.115, 43.077, 48.305, 44.784,
44.934, 46.253, 46.203, 48.36, 47.36, 48.872, 44.803)), .Names = c("Site",
"Day", "Month", "Quality"), class = "data.frame", row.names = c(NA,
-496L))
模型我是这样写的
library(lsmeans)
fit1<-aov(Quality~Site*Month + (Site*Month)/Day, data=test)
那个模型似乎适合我。我了解如何评估站点和月份的交互项和主要影响，但由于某种原因我正在努力评估日期。我试过了
dayeffects<-lsmeans(fit1, pairwise~Site*Month/Day, adjust="bonferroni")
results <- dayeffects[[2]]
summary(results)[!is.na(summary(results)[,4]),]
但这似乎是在测试每个成对比较，而不是遵循嵌套结构。就像我上面说的，我只想比较同一个月和同一站点内发生的天数。
虽然我知道我可以从上面进行我想要的比较，但我觉得我可能做错了什么。它还使 bonferroni 调整过度正确。
任何帮助都会很棒。谢谢
好了，我理清思路了。您可以通过 lsmeans(model, pairwise ~ Day | Site * Month ) 在 Site 和 Month 中比较它们。 （我使用 ~ (Site*Month)/Day 作为模型，但在本主题中，~ Site + Month + Site:Month:Day 或 ~ Site * Month * Day return 结果相同，因为你的 Day 是独一无二的）。
fit <-aov(Quality ~ (Site*Month)/Day, data=test) # this model is equivalent to OP's one.
res <- lsmeans(fit, pairwise ~ Day | Site * Month, adjust="bonferroni")
results <- summary(res[[2]])[!is.na(summary(res[[2]])[,4]),]
> results[25:30,]
Site = H, Month = September:
contrast estimate SE df t.ratio p.value
H6 - H7 -1.0626904 0.5348000 470 -1.987 0.2850
H6 - H8 -0.1373969 0.6588578 470 -0.209 1.0000
H6 - H9 0.3862017 0.5567090 470 0.694 1.0000
H7 - H8 0.9252934 0.6504561 470 1.423 0.9332
H7 - H9 1.4488921 0.5467399 470 2.650 0.0499
H8 - H9 0.5235987 0.6685859 470 0.783 1.0000
检查 p.adjustment
res0 <- lsmeans(fit, pairwise ~ Day | Site * Month, adjust="none")
results0 <- summary(res0[[2]])[!is.na(summary(res0[[2]])[,4]),] # get non-adjusted p.value
res0.p <- p.adjust(results0$p.value[25:30], "bonferroni") # semi-manually p.adjustment
identical(results$p.value[25:30], res0.p) # [1] TRUE
我有一个模型，我正在尝试研究嵌套效应的某些成对比较。我不确定我是否正确编写了模型，而且我不了解如何实际评估嵌套项。
我的数据框有一个名为 'quality' 的响应变量和三个名为 "site" "month" 和 "day" 的预测变量。在我的实验设置中，我测量了每个人的质量。有两个站点。我在 4 个月内对每个网站进行了抽样。每个月有连续 4 天的抽样。我想知道某一天的个体是否与同一个月内其他日子的个体具有显着不同的质量。我对比较一个月的天数和另一个月的天数不感兴趣。
我的数据框如下
test<- structure(list(Site = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("H", "W"), class = "factor"),
Day = structure(c(19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 15L, 15L, 15L,
15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L,
15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 26L, 26L, 26L, 26L,
26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L,
26L, 26L, 8L, 8L, 8L, 8L, 8L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L,
17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L,
17L, 17L, 17L, 17L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L,
14L, 14L, 14L, 14L, 14L, 25L, 25L, 25L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 22L,
22L, 7L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,
16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,
16L, 16L, 16L, 16L, 16L, 16L, 13L, 13L, 13L, 13L, 13L, 13L,
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 24L, 24L, 24L,
24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L,
21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L, 21L,
21L, 21L, 21L, 21L, 21L, 21L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L,
23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L,
23L, 23L, 23L, 23L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L,
20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L,
20L, 20L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L), .Label = c("H1", "H10", "H11", "H12", "H13", "H14",
"H15", "H16", "H2", "H3", "H4", "H6", "H7", "H8", "H9", "W10",
"W11", "W12", "W13", "W2", "W3", "W4", "W6", "W7", "W8",
"W9"), class = "factor"), Month = structure(c(3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 2L, 2L, 2L, 2L, 2L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("August",
"November", "October", "September"), class = "factor"), Quality = c(42.535,
46.651, 45.466, 43.483, 44.896, 46.581, 47.494, 47.529, 46.562,
45.111, 45.982, 48.367, 47.39, 45.388, 46.313, 44.732, 48.641,
46.614, 45.234, 45.96, 44.795, 44.333, 46.559, 46.826, 44.166,
45.19, 46.661, 45.481, 46.828, 43.487, 49.505, 48.558, 45.218,
44.802, 43.975, 47.23, 44.85, 46.213, 44.726, 43.036, 47.211,
45.536, 44.62, 44.297, 36.115, 39.314, 42.349, 44.919, 46.296,
46.317, 45.858, 45.036, 45.861, 48.85, 45.337, 45.03, 47.4,
48.78, 49.829, 45.12, 45.599, 43.235, 44.735, 44.889, 45.666,
46.475, 44.888, 46.215, 42.242, 46.341, 45.992, 43.549, 46.612,
44.232, 42.706, 42.064, 43.837, 43.351, 41.064, 44.364, 42.597,
45.561, 44.51, 45.184, 44.896, 45.772, 47.43, 44.08, 44.697,
45.141, 43.776, 47.175, 46.115, 43.39, 47.426, 47.636, 43.672,
45.987, 45.338, 46.644, 42.192, 47.011, 45.856, 44.764, 36.285,
33.741, 34.324, 35.101, 46.844, 42.52, 48.649, 44.364, 44.688,
45.822, 44.945, 44.311, 44.684, 42.787, 45.516, 46.16, 46.289,
45.661, 45.772, 43.845, 48.717, 46.567, 44.719, 46.585, 45.33,
45.995, 48.053, 44.734, 51.233, 44.597, 45.742, 46.567, 46.478,
44.382, 47.316, 46.205, 45.111, 47.575, 46.014, 44.533, 45.347,
45.983, 47.053, 44.855, 48.021, 45.155, 49.248, 45.634, 48.815,
45.413, 43.091, 47.854, 45.19, 47.495, 47.323, 48.076, 44.183,
43.182, 46.267, 41.58, 44.237, 45.607, 48.517, 44.639, 44.773,
42.787, 43.965, 46.629, 46.256, 47.688, 44.126, 44.712, 47.097,
44.561, 47.306, 45.323, 46.328, 45.832, 46.075, 46.778, 47.445,
45.582, 47.691, 45.193, 48.453, 46.301, 44.847, 43.675, 46.066,
47.896, 45.2, 44.959, 47.401, 46.267, 45.743, 47.411, 46.926,
46.24, 46.212, 44.988, 36.552, 38.027, 47.355, 40.147, 38.094,
39.043, 37.589, 46.491, 46.413, 43.92, 45.228, 46.319, 44.764,
47.376, 43.924, 45.203, 45.418, 45.684, 46.34, 43.655, 44.365,
46.927, 48.269, 45.473, 46.451, 42.752, 48.346, 47.832, 46.534,
46.47, 43.282, 47.749, 44.856, 46.551, 45.925, 45.669, 47.263,
44.367, 47.017, 42.922, 44.904, 48.85, 45.535, 48.512, 46.154,
47.306, 46.571, 46.619, 46.092, 43.808, 47.7, 48.482, 44.407,
45.442, 44.771, 46.373, 47.777, 43.012, 46.154, 45.203, 46.443,
43.461, 45.714, 40.776, 48.949, 45.72, 48.269, 45.782, 43.945,
45.382, 43.729, 44.187, 45.267, 46.012, 42.234, 43.431, 41.973,
45.597, 45.993, 46.303, 44.493, 44.981, 46.487, 45.01, 47.009,
46.904, 48.277, 48.585, 48.625, 47.511, 44.011, 42.21, 47.124,
44.244, 47.76, 47.299, 45.278, 45.564, 44.621, 46.75, 45.396,
44.947, 46.185, 45.399, 46.095, 49.545, 47.211, 43.613, 48.494,
44.102, 45.888, 45.473, 47.222, 46.681, 45.863, 47.834, 48.386,
46.979, 46.318, 46.061, 46.347, 47.976, 47.079, 48.254, 47.643,
46.244, 46.717, 44.574, 45.177, 44.879, 46.485, 47.416, 50.235,
45.626, 48.117, 44.529, 44.281, 47.087, 47.356, 43.234, 45.841,
43.487, 42.997, 35.322, 45.554, 44.973, 43.396, 43.023, 44.65,
47.088, 41.934, 45.704, 44.559, 37.969, 42.687, 42.995, 45.287,
45.21, 43.335, 46.892, 45.534, 44.19, 43.606, 44.173, 49.334,
44.888, 47.477, 47.054, 41.041, 46.629, 45.049, 44.478, 40.278,
43.044, 43.575, 46.194, 42.688, 41.361, 46.828, 45.534, 47.395,
45.431, 45.433, 45.331, 43.947, 47.371, 48.308, 45.726, 41.833,
45.782, 44.756, 45.406, 45.661, 43.447, 46.932, 45.495, 44.349,
40.493, 43.603, 48.151, 44.037, 44.379, 45.934, 44.854, 42.321,
46.198, 44.622, 46.077, 45.306, 48.951, 47.972, 42.581, 43.608,
45.988, 44.955, 45.097, 46.768, 44.722, 45.971, 46.612, 48.956,
47.669, 47.757, 47.189, 44.184, 48.464, 49.546, 48.021, 45.448,
45.573, 46.778, 45.769, 45.419, 45.277, 47.489, 46.762, 46.238,
47.509, 47.249, 46.243, 46.124, 46.801, 47.385, 43.614, 44.661,
45.96, 48.791, 47.872, 42.402, 45.651, 45.927, 43.781, 49.923,
47.153, 46.87, 43.767, 47.3, 46.897, 44.932, 45.135, 50.124,
45.366, 45.063, 45.958, 46.731, 43.863, 45.095, 47.755, 45.446,
45.145, 45.998, 46.377, 44.369, 46.485, 48.852, 45.365, 45.934,
44.856, 48.195, 45.424, 49.05, 46.115, 43.077, 48.305, 44.784,
44.934, 46.253, 46.203, 48.36, 47.36, 48.872, 44.803)), .Names = c("Site",
"Day", "Month", "Quality"), class = "data.frame", row.names = c(NA,
-496L))
模型我是这样写的
library(lsmeans)
fit1<-aov(Quality~Site*Month + (Site*Month)/Day, data=test)
那个模型似乎适合我。我了解如何评估站点和月份的交互项和主要影响，但由于某种原因我正在努力评估日期。我试过了
dayeffects<-lsmeans(fit1, pairwise~Site*Month/Day, adjust="bonferroni")
results <- dayeffects[[2]]
summary(results)[!is.na(summary(results)[,4]),]
但这似乎是在测试每个成对比较，而不是遵循嵌套结构。就像我上面说的，我只想比较同一个月和同一站点内发生的天数。
虽然我知道我可以从上面进行我想要的比较，但我觉得我可能做错了什么。它还使 bonferroni 调整过度正确。
任何帮助都会很棒。谢谢
好了，我理清思路了。您可以通过 lsmeans(model, pairwise ~ Day | Site * Month ) 在 Site 和 Month 中比较它们。 （我使用 ~ (Site*Month)/Day 作为模型，但在本主题中，~ Site + Month + Site:Month:Day 或 ~ Site * Month * Day return 结果相同，因为你的 Day 是独一无二的）。
fit <-aov(Quality ~ (Site*Month)/Day, data=test) # this model is equivalent to OP's one.
res <- lsmeans(fit, pairwise ~ Day | Site * Month, adjust="bonferroni")
results <- summary(res[[2]])[!is.na(summary(res[[2]])[,4]),]
> results[25:30,]
Site = H, Month = September:
contrast estimate SE df t.ratio p.value
H6 - H7 -1.0626904 0.5348000 470 -1.987 0.2850
H6 - H8 -0.1373969 0.6588578 470 -0.209 1.0000
H6 - H9 0.3862017 0.5567090 470 0.694 1.0000
H7 - H8 0.9252934 0.6504561 470 1.423 0.9332
H7 - H9 1.4488921 0.5467399 470 2.650 0.0499
H8 - H9 0.5235987 0.6685859 470 0.783 1.0000
res0 <- lsmeans(fit, pairwise ~ Day | Site * Month, adjust="none")
results0 <- summary(res0[[2]])[!is.na(summary(res0[[2]])[,4]),] # get non-adjusted p.value
res0.p <- p.adjust(results0$p.value[25:30], "bonferroni") # semi-manually p.adjustment
identical(results$p.value[25:30], res0.p) # [1] TRUE