Evaluating pairwise comparisons of nested effects

Question

I have a model which I'm attempting to investigate some pairwise comparisons of nested effects. I'm not sure if I've written the model correctly, and I'm not understanding how to actually evaluate the nested term.

My data frame has one response variable called 'quality' and three predictor variables called "site" "month" and "day". In my experimental setup I measured quality on each individual. There were two sites. I sampled each site over 4 months. Each month had 4 consecutive days of sampling. I would like to know if individuals on one day are of significantly different quality to individuals from other days within the same month. I'm not interested in comparing days to one month to days from another month.

My data frame is as follows

   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))

I've written the model like this

    library(lsmeans)
fit1<-aov(Quality~Site*Month + (Site*Month)/Day, data=test)

That model seems to work for me. I understand how to evaluate the interaction term and main effects of site and month, but I'm struggling to evaluate day for some reason. I've tried

dayeffects<-lsmeans(fit1, pairwise~Site*Month/Day, adjust="bonferroni")
results <- dayeffects[[2]]
summary(results)[!is.na(summary(results)[,4]),]

But this appears to test every pairwise comparison, rather than following the nesting structure. Like I said above, I only want to compare days that occur within the same month and site.

While I know that I could just take the comparisons I want from above, I feel like I'm probably doing something wrong. Also it makes the bonferroni adjustment overcorrect.

Any help would be great. Thanks

Sorry, I probably should add that I do want to evaluate the other effects in the model, hence why I've written the model that way. I just don't have any issues that part of the analysis — Conor Neilson, Oct 09 '16 at 23:26
Yep but I want to test each site specifically, not average over that variable. I coded it that way because I thought you needed unique IDs in nested designs? — Conor Neilson, Oct 10 '16 at 00:00
@cuttlefish44 Thanks! I used `test<-lsmeans(fit2, pairwise~Month/Day, adjust="bonferroni")` `results <- test[[2]]` `summary(results)[!is.na(summary(results)[,4]),]` And that shows only the tests of interest. But it still reports that the bonferrini is being adjusted for 5356 tests. Do you know why?/how to do the correct adjustment? — Conor Neilson, Oct 10 '16 at 00:25

score 0 · Accepted Answer · answered Oct 10 '16 at 04:03

OK, I cleared my thoughts. You can compare them within Site and Month by lsmeans(model, pairwise ~ Day | Site * Month ). (I used ~ (Site*Month)/Day as model, but in this topic, ~ Site + Month + Site:Month:Day or ~ Site * Month * Day return the same results because your Day is unique).

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

check 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

Evaluating pairwise comparisons of nested effects

1 Answers1