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

• Having even a lot of 0 values doesn't necessarily make a variable zero-inflated, so I wouldn't worry about that for now. More important for this study: do the times for the 3 behaviors always add up to 60 for each row/observation, as they do for these examples? I think that there's a reasonably simple way to model these data, especially if that's the case. Please provide that information by editing the question itself, as comments are easy to overlook and can be deleted.
– EdM
Commented Nov 1, 2022 at 14:19
• Hi, yes the behaviours each hour always add to 60, and time goes from 0 to 23. and the distribution of the data definitelly shows lots of 0 for each behaviour. 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 alocated in each paddock (ns and sh), then at 2 days we swaped the groups. A Commented Nov 2, 2022 at 5:04
– EdM
Commented Nov 7, 2022 at 10:35
• If you actually have data on the cow's status at baseline, I would recommend just fitting a transition state model - i.e. estimate the Markov Chain of what-my-cow-is-doing. Commented Nov 15, 2022 at 16:15
• The answer with simulated data now includes a suggestion for how to model the transitions between behaviors instead of the behavior probabilities, if you have minute-by-minute data.
– EdM
Commented Nov 25, 2022 at 16:26

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

• I was thinking about testing cows if they differ between them in each treatment and then use cow as a repetition, instead of repeated measurement, then run a principal component analysis. The questions that I want to ask are: Are the hours affecting the behaviours on the treatments? Also, is s and ns different in relation to the beviours? I manage to make the zero inflated model to work, but didnt put the hours as factor, just as a fixed variable. So the output wasnt what I needed. I am cracking my head to analyse it. Commented Nov 3, 2022 at 5:04
• @jacqueskrtickaCarvalho I think that my suggestion for multinomial or ordinal regression makes the most sense for answering the question you pose of the data, but those can be closely related to Poisson-type models. It's hard to know what the problems with your models are unless you post the code that you used. Please edit the original question to add the code that you tried and the results that you got.
– EdM
Commented Nov 4, 2022 at 16:19
• There was no problem with the model, it just didnt answer my questions. Is just there are differences per hour in each behaviour comparing s and ns? (ex. at 13 pm cattle significantly spent more time active on shelterbelt. Commented Nov 7, 2022 at 6:15
• So I tried to use the gee model but I had this error: fitmod <- ordLORgee(formula =Active ~ factor(time) + factor(trt) + factor(day), data = datam, id = cow, LORstr = "uniform", repeated = time) It results in: Error in ordLORgee(formula = datam $Active ~ factor(time) + factor(trt) + : 'repeated' does not have unique values per 'id' What am I doing wrong? Commented Nov 8, 2022 at 3:26 • @jacqueskrtickaCarvalho I don't think that you need the repeated argument here; it won't work with your time values anyway, as you seem to have 4 sets of the 24 time values for each cow. If you want/need to use repeated, it needs to be a series of unique number for a cow, in your case from 1 to 96 if you have data for 4 days. The first observation would be 1, the last, 96, for 4 days with 24 observations per day. You must be very careful in setting up the data; I don't see that your formula provides data on the behaviors other than Active. Check the help page and the examples. – EdM Commented Nov 8, 2022 at 13:35 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.