I have zero inflated data, with discrete variables. Is it possible to use zero inflated poisson model? I have behaviour data to analyse: the time budgets of cows (three columns of 3 different behaviours per hour). Values range from 0 to 60 minutes each hour (e.g. 30 min active, 10 min lying and 20 standing). The behaviours each hour always add to 60, and time goes from 0 to 23.
For variables I also have day and cows (which goes from 1 to 20, and I add on the model as a random variable). it is a crossover design, meaning that all experimental units were exposed to both treatments (shelterbelts and non-shelterbelts). So, 20 cows in two groups of ten were allocated in each paddock (ns and sh), then at 2 days we swapped the groups.
The distribution of the data is zero inflated for each behaviour, so can I use the zero inflated poisson model to find differences between treatments and test interaction between time and treatment in relation to behaviour? note:Data ranges from 0 to 60, so can I use this model for discrete data?
ex:
trt Active Lying Standing  time
ns   10      20     30       0
s     30     20      10      1
s     50     10      0       2
ns    0      50      10      3

I want to find if there is difference between shelterbelt and non shelterbelt generally in relation to the time budgets, also if there is an effect of hour of the day on the behaviour in each treatment.

As my other answer is already long, I'll use the actual data linked in a comment to illustrate circular splines for handling hourly data and further analysis of the time by treatment interaction. Otherwise, the approach is as in my other answer for these data aggregated by the hour.
I converted the linked data into a data frame called datam; the format is similar to that shown in the question but includes an indentifier for cow. Several packages help here:
library(tidyr) # formatting data
library(nnet)  # multinomial fit
library(pbs)   # for circular spline
library(aod)   # for Wald test
library(emmeans)
library(ggplot2)

Data were converted into two other forms: first, a long form that puts each Behavior on a separate line for each hour, with the corresponding "count" values placed in a minutes column; then, a nested format that keeps all data for a cow in a single row, to simplify bootstrapping.
## put datam in long form, separating behaviors
## simplifies use of other packages
datamLong <- pivot_longer(datam,
  cols=c("Active","Lying",Standing), 
  names_to="Behavior", values_to="minutes")
## also nested by cow for cluster bootstrap sampling
dfNest <- nest(datamLong,data=!cow)

The pbs package simplifies use of circular splines to fit the multinomial model to the data. That package was recommended in an answer by @kjetilbhalvorsen to another question. The splines are restricted so that the values at 0 and 24 hours are the same.
mn1 <- multinom(Behavior ~
  trt * pbs(time,df=4,Boundary.knots=c(0,24)),
  data=datamLong, weights=minutes, trace=FALSE)

The two groups in the crossover design (sh -> ns; ns -> sh; as identified from the data) are sampled separately in the bootstrapping by cow to get the robust covariance matrix.
B = 200 # number of samples
b = 0 # for successful samples
# p is number of parameters
# put into the order used by multinom
p = length(c(t(coef(mn1)))) 
bar <- rep(0, p)
coefSet <- NULL # for keeping individual estimates
cov <- matrix(0, nrow = p, ncol = p) # for covariance estimate

set.seed(204)
## sample sh -> ns and ns -> sh groups separately
for (boots in 1:B) {
bootSam <- dfNest[match(c(
  sample(c(2,4,5,20,22,23,24,25,34,38),
    size=10,replace=TRUE),
  sample(c(8,19,26,31,32,35,36,37,39,40),
    size=10,replace=TRUE)),dfNest$cow),]
## unnest for multinomial modeling
bootSam <- unnest(bootSam,col=!cow)
bootMult <- tryCatch(multinom(Behavior ~ 
  trt * pbs(time,df=4,Boundary.knots=c(0,24)),
  data=bootSam, weights=minutes, trace=FALSE),
  error = function(...) list(fail = TRUE))
if (length(bootMult$fail) && bootMult$fail) 
                next
b = b + 1
cof <- c(t(coef(bootMult)))
coefSet <- cbind(coefSet,c(t(coef(bootMult))))
bar <- bar + cof
cof <- as.matrix(cof)
cov <- cov + cof %*% t(cof)
}
bar <- bar/b
boot.coef <- bar
bar <- as.matrix(bar)
cov <- (cov - b * bar %*% t(bar))/(b - 1L)

The Wald test (from the aod package) based on the entire set of interaction coefficients and the robust covariance matrix cov addresses whether there is an interaction between treatment and time with respect to outcome. The Terms argument represents the numerical order of the interaction terms among all the coefficients in c(t(coef(mn1))), which is also their order in the covariance matrix. If you fit a model with a different number of parameters to estimate, you will have to change that specification.
# uses aod package Wald test on the interaction coefficients
# set Terms to represent the  interaction coefficient positions
wald.test(cov,c(t(coef(mn1))),Terms=c(7:10,17:20))
# Wald test:
# ----------
# 
# Chi-squared test:
# X2 = 114.7, df = 8, P(> X2) = 0.0

There is a statistically significant interaction.
To evaluate the practical importance of the statistically significant interaction, next use the emmeans package to determine values of the model estimates and of the differences between sh and ns over time for each Behavior. The vcov. argument specifies use of the robust covariance matrix cov.
rg1 <- ref_grid(mn1,vcov.=cov,at=list(time=seq(0,24,by=0.25)))
emm1 <- emmeans(rg1,~trt.vs.ctrl1~trt|time+Behavior)
df1 <- data.frame(emm1$emmeans)
df2 <- data.frame(emm1$contrast)

There might be a more elegant way to do this in emmeans, but this method works.
To plot behavior probabilities over time, treatments distinguished by color, with 95% CI based on the robust covariance matrix:
ggplot(df1,mapping=aes(x=time,y=prob,group=trt)) +
  geom_line(aes(color=trt)) +
  geom_ribbon(aes(ymin=lower.CL,ymax=upper.CL),
  alpha = 0.1) +
  facet_wrap(facets=vars(Behavior)) + 
  scale_x_continuous(breaks=seq(0,24,by=6))


The df2 data frame derived from emm1$contrast includes sh-ns estimates and p-values for all behaviors and times at 0.25-hour intervals. You can interrogate that as illustrated in the other answer. To plot treatment differences over time for each behavior, again with 95% CI:
ggplot(df2,mapping=aes(x=time,y=estimate)) + 
  geom_line() +facet_wrap(facets=vars(Behavior)) + 
  geom_ribbon(aes(ymin=estimate-1.96*SE,
    ymax=estimate+1.96*SE),alpha=0.1) +
  ggtitle("sh minus ns") +
  ylab("Probability difference") + 
  scale_x_continuous(breaks=seq(0,24,by=6))


The treatment differences are typically within $\pm 0.1$ on the probability scale. Apply your understanding of the subject matter to determine the practical significance.
