Outcome to model
A cow has 3 mutually exclusive choices at any time: lying, standing, and active. You record the choice over time as the outcome.
As AdamO suggests in a comment, this might most efficiently be evaluated with a transition-state model of factors that affect a cow's decision to move from one behavior to another. That, however, would require having the minute-by-minute observations available for all animals, annotated by time of day and treatment. If you have such data, then you could use the multi-state modeling tools in the R msm package to model the probabilities of transitions between states as a function of treatment and time of day. That could be more informative than simply modeling the probability of each state over time.
You show data already aggregated by the hour instead. The most natural way to model such data is a multinomial or ordinal regression model, extensions of binomial regression to multiple categorical outcomes. The choice between multinomial and ordinal models is based on your understanding of whether there is a natural ordering of the outcomes. The choice between either of those and a model of transition probabilities (if you have the necessary data) depends on how you want to use the model and interpret the effects of time-of-day and treatment.
A Poisson model has no upper limit to the counts, while your values can't exceed 60. There is a so-called "Poisson trick" that can sometimes be used to model multinomial responses. That doesn't, however, seem to be what you used. It requires a special data setup, and it doesn't extend well to correlated data like yours.
I'll focus on multinomial modeling. I set up a data set that matches the structure of your data; code is at the end. It was simplest to reshape your data format into separate rows for each behavior at each combination of predictors, with the associated number of counts (from 0 to 60) noted. That makes behavior the outcome, with observations weighted by the corresponding counts.
The multinom() function in the R nnet package posed the fewest problems of methods I tried. That uses one behavior as the reference (I chose lying) and models the log-odds of each of the other behaviors relative to that.
Predictors in the model
The regression model must include treatment as the predictor of main interest; I called it tx. To see if time-of-day (tod) alters the effect of tx, include an interaction between tx and tod.
Handling of time requires some thought. Treating time as a categorical predictor requires estimating many coefficients.
library(nnet)
nn1 <- multinom(behavior ~ tx * factor(tod),
data=dfCowLong, weights=counts, trace=FALSE)
dim(coef(nn1))
# [1] 2 48
You have 48 coefficients for the log-odds of each of the two behaviors standing and active, versus the baseline of lying
Treating time as a circular predictor could be better, as explained on this page. You model time (appropriately scaled) as a sum of a sine and a cosine term, which cuts down the number of coefficients dramatically while still allowing for an interaction with trt.
nn2 <- multinom(behavior ~ tx *
(sin(2 * pi * tod/nTimePerDay) +
cos(2 * pi * tod/nTimePerDay)),
data=dfCowLong, weights=counts, trace=FALSE)
## scales by nTimePerDay from data setup
dim(coef(nn2))
## [1] 2 6
That cuts the number of coefficients by a factor of 8. A regression spline for time might also be used.
Those models, however, assume that the observations are independent. They aren't, as the same cows were monitored repeatedly.
Handling repeated measurements on cows
Bootstrapping can handle the correlations of within-cow observations by providing a type of robust estimate of the coefficient covariance matrix. This is similar to the approach of generalized estimating equations, as it makes a working assumption of independence to get the coefficient estimates but then adjusts the covariance matrix.
After fitting your model on the full data set to get the coefficient estimates, you re-fit it on multiple bootstrapped samples of your data, resampling by cow ("cluster bootstrap") rather than by individual observation.
The following uses the nest() function in the tidyr package to collect the data by id. Bootstrap by id (separately for the two halves of the cows, to respect the cross-over structure), unnest(), then collect data from each bootstrap sample. This borrows code snippets from Frank Harrell's bootcov() function in the rms package. This takes a while to run for 200 bootstraps; you could alter the number provided by B in the code If you're impatient.
B = 200 ## number of bootstraps
## initialize variables for collecting
b = 0 ## for successful bootstraps
## get numbers of coefficients
## multinom coefficient are a matrix
p = length(c(t(coef(nn1))))
bar <- rep(0, p) ## for average coifs
cov <- matrix(0, nrow = p, ncol = p) ## for coefficient cross-products
## nest to collect by id
dfNest <- nest(dfCowLong,data=!id)
halfSample <- as.integer(nCows/2) ## for 2 groups of cows
set.seed(204)
## do the sampling and collect results
for (boots in 1:B) {
bootSam <- dfNest[c(sample(1:halfSample,size=halfSample,replace=TRUE),sample(halfSample+1:nCows,size=halfSample,replace=TRUE)),]
bootSam <- unnest(bootSam,col=!id)
bootMult <- tryCatch(multinom(behavior ~ tx * factor(tod), weights=counts, trace=FALSE, data=bootSam), error = function(...) list(fail = TRUE))
## jump to next sample if problem with fit
if (length(bootMult$fail) && bootMult$fail)
next
## update if successful
b = b + 1
cof <- c(t(coef(bootMult)))
bar <- bar + cof
cof <- as.matrix(cof)
#outer product of coefficient estimates
cov <- cov + cof %*% t(cof)
}
## after all bootstraps
bar <- bar/b ## average coefficients over successful bootstraps
boot.coef <- bar ## save separately
bar <- as.matrix(bar)
## following is analogous to the formula for a variance of a scalar
cov <- (cov - b * bar %*% t(bar))/(b - 1L)
Then you can use the tools of the emmeans package to inquire about the results for any combinations of predictor values. It accepts a vcov. argument that replaces the model's original coefficient covariance matrix, to account for the lack of independence. It reports predictions directly in terms of probability. Here's an example:
library(emmeans)
emm1 <- emmeans(nn1,trt.vs.ctrl1~tx|behavior+tod, vcov.=cov)
pred.df1 <- data.frame(summary(emm1)$emmeans)
A plot of modeled behavior probabilities over time, as functions of treatment, with 95% confidence limits:
library(ggplot2)
ggplot(data=pred.df1,mapping=aes(x=tod,y=prob,group=tx)) +
geom_line(aes(color=tx)) +
scale_x_continuous(breaks=seq(0,24,by=6)) +
xlab("Time of day") + ylab("Behavior probability") +
facet_wrap(facets=vars(behavior)) +
geom_ribbon(mapping=aes(ymin=lower.CL,ymax=upper.CL),alpha=0.1)
With emmeans you can compare estimates for any combinations of time and treatment. The emm1 object produced above includes a contrast object that provides the treatment differences and p-values for all combinations of behaviors and times. It might be easiest to convert that to a data frame, which you then can interrogate. For example, to get the difference in treatment effects for all behaviors at tod=23:
pred.df2 <- data.frame(emm1$contrast)
names(pred.df2)
# [1] "contrast" "behavior" "tod" "estimate" "SE" "df" "t.ratio" "p.value"
print(pred.df2[pred.df2$tod==23,],digits=3)
# contrast behavior tod estimate SE df t.ratio p.value
# 67 1 - 0 lying 23 -0.0779 0.0132 96 -5.90 5.46e-08
# 68 1 - 0 standing 23 0.0434 0.0144 96 3.02 3.24e-03
# 69 1 - 0 active 23 0.0345 0.0145 96 2.38 1.95e-02
Those p-values aren't corrected for multiple comparisons among behaviors and treatments, however. A conservative (but simple) Bonferroni correction would limit "significance" at a family-wise error rate of 0.05 to p-values less than 0.00069 = 0.05/(24*3).
Other options
As an alternate type of robust covariance estimate, you might consider a marginal model that uses the generalized estimating equation approach to handle the within-cow correlations. The R multgee package is designed for that; it handles both multinomial and ordinal outcomes. It seems to expect individual observations of outcomes over time rather than aggregated outcome data. If you have minute-by-minute observations for each animal, that could be a reasonable choice for a multinomial model (although the package is very strict about how it expects the data to be arranged).
The R mgcv package can handle multinomial outcomes, cyclic regression splines to handle time flexibly, and random effects. I didn't readily get multgee or mgcv to work, however, with aggregated data. In contrast, bootstrapping with multinom() worked well. The R mixed models task view has additional suggestions for modeling correlated multinomial or ordinal outcomes.
Code to construct the data set.
nCows <- 20 ## should be even
nTimePerDay <- 24
nDays <- 4 ## should be even
nObs <- nTimePerDay * nDays
## organize by cow
id <- rep(1:nCows,each=nObs)
## first half of cows get one tx, rest the other
## switch half-way through
tx=c(rep(c(rep(0,nObs/2),rep(1,nObs/2)),nCows/2),rep(c(rep(1,nObs/2),rep(0,nObs/2)),nCows/2))
## set up times of day to match the above
tod <- rep(1:nTimePerDay,nCows * nDays)
# put together
dfCow <- data.frame(id=id,tx=tx,tod=tod)
## random baseline lying-down probs
set.seed(203)
idBase <- rnorm(nCows,mean=0.5,sd=0.1)
## define probability differences over time:
addT <- c(seq(0,0.1,length.out=nTimePerDay/2),rev(seq(0,0.1,length.out=nTimePerDay/2)))
## sum up baseline, subtract 0.1 for tx=1, add time-of-day effect
dfCow[,"lyingP"] <- idBase[dfCow$id] - 0.1*dfCow$tx + addT[dfCow$tod]
## convert tx to factor
dfCow$tx <- factor(dfCow$tx)
## function to make 60 multinomial samples based on lyingP
## use values for lying down, then equal probs for other2 choices
do60 <- function(x) rmultinom(1,size=60,prob=c(x,(1-x)/2,(1-x)/2))
do60 <- Vectorize(do60)
countVals<- t(do60(dfCow[,"lyingP"]))
colnames(countVals) <- c("lying","standing","active")
## append counts (each summing to 60) to data
dfCow <- cbind(dfCow,countVals)
## long-form data play better with other packages
library(tidyr)
dfCowLong <- pivot_longer(dfCow,cols=c("lying","standing","active"),names_to="behavior",values_to="counts")
dfCowLong$behavior <- factor(dfCowLong$behavior, levels=c("lying","standing","active"))
