Évaluation des comparaisons par paires d'effets imbriqués

Question

J'ai un modèle dans lequel j'essaie d'étudier des comparaisons par paires d'effets imbriqués. Je ne suis pas sûr d'avoir correctement écrit le modèle et je ne comprends pas comment évaluer le terme imbriqué.

Ma trame de données a une variable de réponse appelée «qualité» et trois variables de prédicteur appelées «site» «mois» et «jour». Dans ma configuration expérimentale, j'ai mesuré la qualité de chaque individu. Il y avait deux sites. J'ai échantillonné chaque site pendant 4 mois. Chaque mois a eu 4 jours consécutifs d'échantillonnage. J'aimerais savoir si les personnes d'un jour sont de qualité sensiblement différente de celles des autres jours du même mois. Je ne suis pas intéressé à comparer des jours à un mois à des jours d'un autre mois.

Mon trame de données est la suivante

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

J'ai écrit le modèle comme celui-ci

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

Ce modèle semble fonctionner pour moi. Je comprends comment évaluer le terme d'interaction et les principaux effets du site et du mois, mais j'ai du mal à évaluer le jour pour une raison quelconque. J'ai essayé

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

Mais cela semble tester chaque comparaison par paire plutôt que de suivre la structure d'imbrication. Comme je l'ai dit ci-dessus, je veux seulement comparer les jours qui se produisent dans le même mois et le même site.

Même si je sais que je pourrais prendre les comparaisons que je veux d'en haut, j'ai l'impression de faire quelque chose de mal. En outre, il rend l'ajustement bonferroni excessif.

Toute aide serait géniale. Merci

Answer 1

OK, j'ai éclairci mes pensées. Vous pouvez les comparer entre Site et Month par lsmeans(model, pairwise ~ Day | Site * Month). (J'ai utilisé ~ (Site*Month)/Day comme modèle, mais dans ce sujet, ~ Site + Month + Site:Month:Day ou ~ Site * Month * Day renvoyer les mêmes résultats parce que votre Day est 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