Repeated measures lmer with groups

Question

Background: I have data that I'd like to fit a lmer model to and a snippet of the data is available below. Basically, predictions were made on the grazing duration Graz_time of cattle in 3 groups of 7 cows. Each group was monitored for 3 days before moving on to the next group for another 3 days etc.

data <- structure(list(Cow = structure(c(1L, 3L, 5L, 7L, 8L, 9L, 12L, 
1L, 3L, 5L, 7L, 8L, 9L, 12L, 1L, 3L, 5L, 7L, 8L, 9L, 12L, 10L, 
14L, 15L, 16L, 18L, 20L, 21L, 10L, 14L, 15L, 16L, 18L, 20L, 21L, 
10L, 14L, 15L, 16L, 18L, 20L, 21L, 2L, 4L, 6L, 11L, 13L, 17L, 
19L, 2L, 4L, 6L, 11L, 13L, 17L, 19L, 2L, 4L, 6L, 11L, 13L, 17L, 
19L), .Label = c("5", "15", "31", "55", "68", "74", "78", "84", 
"115", "162", "163", "197", "266", "271", "272", "288", "292", 
"391", "430", "449", "756"), class = "factor"), Group = structure(c(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, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("1", 
"2", "3"), class = "factor"), Day = structure(c(1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("1", 
"2", "3"), class = "factor"), BCS = structure(c(6L, 6L, 4L, 3L, 
5L, 3L, 5L, 6L, 6L, 4L, 3L, 5L, 3L, 5L, 6L, 6L, 4L, 3L, 5L, 3L, 
5L, 4L, 6L, 6L, NA, 6L, 1L, 3L, 4L, 6L, 6L, NA, 6L, 1L, 3L, 4L, 
6L, 6L, NA, 6L, 1L, 3L, 2L, 5L, 5L, 6L, 6L, 6L, 4L, 2L, 5L, 5L, 
6L, 6L, 6L, 4L, 2L, 5L, 5L, 6L, 6L, 6L, 4L), .Label = c("1.5", 
"1.75", "2", "2.25", "2.5", "2.75"), class = "factor"), Parity = 
structure(c(2L, 
1L, 2L, 2L, 2L, 3L, 3L, 2L, 1L, 2L, 2L, 2L, 3L, 3L, 2L, 1L, 2L, 
2L, 2L, 3L, 3L, 5L, 1L, 1L, 1L, 1L, 1L, 3L, 5L, 1L, 1L, 1L, 1L, 
1L, 3L, 5L, 1L, 1L, 1L, 1L, 1L, 3L, 6L, 1L, 3L, 2L, 1L, 1L, 4L, 
6L, 1L, 3L, 2L, 1L, 1L, 4L, 6L, 1L, 3L, 2L, 1L, 1L, 4L), .Label = c("1", 
"2", "3", "4", "5", "6"), class = "factor"), Month_lact = structure(c(1L, 
4L, 2L, 2L, 2L, 4L, 3L, 1L, 4L, 2L, 2L, 2L, 4L, 3L, 1L, 4L, 2L, 
2L, 2L, 4L, 3L, 7L, 6L, 6L, 7L, 7L, 6L, 7L, 7L, 6L, 6L, 7L, 7L, 
6L, 7L, 7L, 6L, 6L, 7L, 7L, 6L, 7L, 8L, 7L, 7L, 7L, 7L, 5L, 4L, 
8L, 7L, 7L, 7L, 7L, 5L, 4L, 8L, 7L, 7L, 7L, 7L, 5L, 4L), .Label = c("1", 
"2", "3", "4", "6", "7", "8", "9"), class = "factor"), L_NL = 
structure(c(2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 
1L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("L", 
"NL"), class = "factor"), Grass = c(2396L, 2396L, 2396L, 2396L, 
2396L, 2396L, 2396L, 2756L, 2756L, 2756L, 2756L, 2756L, 2756L, 
2756L, 3451L, 3451L, 3451L, 3451L, 3451L, 3451L, 3451L, 2863L, 
2863L, 2863L, 2863L, 2863L, 2863L, 2863L, 2532L, 2532L, 2532L, 
2532L, 2532L, 2532L, 2532L, 2358L, 2358L, 2358L, 2358L, 2358L, 
2358L, 2358L, 3211L, 3211L, 3211L, 3211L, 3211L, 3211L, 3211L, 
2829L, 2829L, 2829L, 2829L, 2829L, 2829L, 2829L, 2552L, 2552L, 
2552L, 2552L, 2552L, 2552L, 2552L), FieldSize_ha = c(3.14, 3.14, 
3.14, 3.14, 3.14, 3.14, 3.14, 1.64, 1.64, 1.64, 1.64, 1.64, 1.64, 
1.64, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 1.59, 1.59, 1.59, 1.59, 
1.59, 1.59, 1.59, 5.92, 5.92, 5.92, 5.92, 5.92, 5.92, 5.92, 10.41, 
10.41, 10.41, 10.41, 10.41, 10.41, 10.41, 2.85, 2.85, 2.85, 2.85, 
2.85, 2.85, 2.85, 2.85, 2.85, 2.85, 2.85, 2.85, 2.85, 2.85, 2.25, 
2.25, 2.25, 2.25, 2.25, 2.25, 2.25), Graz_time = c(315L, 444L, 
273L, 426L, NA, 381L, 486L, 369L, 345L, 276L, 348L, 297L, 363L, 
474L, 354L, 375L, 288L, 312L, 447L, 342L, 444L, 291L, 270L, 303L, 
120L, 189L, 426L, 324L, 285L, 531L, 483L, 447L, 366L, 393L, 525L, 
435L, 483L, 447L, 459L, 417L, 651L, 645L, 480L, 462L, 573L, 486L, 
543L, 270L, 288L, 366L, 207L, 327L, 351L, 372L, 399L, 312L, 531L, 
447L, 225L, 528L, 408L, 330L, 363L)), .Names = c("Cow", "Group", 
"Day", "BCS", "Parity", "Month_lact", "L_NL", "Grass", "FieldSize_ha", 
"Graz_time"), row.names = c(NA, -63L), class = "data.frame")

Having read THIS answer I can see that if I construct the following model mod that the auto-correlation expected within Cow will be accounted for. I also want to account for the variance that may be present as a result of the Group category also and so I think that this also accounts for that, i.e. Cow is nested within Group.

Question: My question is then, do i need to add Day as a fixed effect in the model also with it already accounted for in effect in Cow? Biologically I am not that interested in Day other than that it is a repeated measure design. I have also read THIS post which I think suggests in this instance that I may not need Day as a fixed effect.

Further info: For some reason I now get an error with this data (it is a snippet of what I have) but some advice on the question would be very helpful.

mod <- lmer(Graz_time ~ BCS + Parity + Month_lact + L_NL + Grass + FieldSize_ha + (1|Group/Cow), data = data, REML=FALSE)

Answer 1

I think your intuition about the inclusion (or exclusion) of Day is correct. Given there is a single Cow measurements per day as well as Day not being of clear biological interest, excluding Day is not problematic despite being part of the repeated measurement design.

Having said the above, the problem with this model is that it is over-parametrised for this snippet of data. For starters, even very lenient rules-of-thumb for multiple regression would have trouble justifying 6 explanatory variables with just 63 observations. Here we have 6 explanatory variables plus a nested random effects design; this seems a bit of stretch. To that extent, certain levels of the nominal variables used as fixed effects exhibit perfect correlation (the triplet of variables Month_lact9, Parity6 and BCS1.75) or in other cases, numerical variables seem quite related (the pair Grass and FieldSize_ha has a Spearman rank correlation of ~ -0.65). These findings are not catastrophic but also raise some worries about how much additional information we get from each new variable. As a consequence, somewhat unsurprisingly, any variation due to random effects in terms of Group, Cow or Group/Cow is unfortunately zero when using REML=TRUE. Given the above, I think one should examine how to simplify the modelling assumptions and/or obtain more data.