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