A CASE STUDY ON ANIMAL BEHAVIOR ANALYSIS USING GAMLSS

Animal behavior studies usually produce large amounts of data and a wide variety of data structures, including nonlinear relationships, interaction effects, nonconstant variance, correlated measures, overdispersion, and zero inflation, among others. We aimed to explore here the potential of generalized additive models for location, scale and shape (GAMLSS) in analyzing data from animal behavior studies. Data from 20 Romane ewes from two genetic lineages submitted to brushing by a familiar observer were analyzed. Behavioral responses through ear posture changes, a count random variable, and the proportion of time to perform the horizontal ear posture, a continuous random variable on the interval (0,1), with non-null probabilities in zero and one, were analyzed. The Poisson, negative binomial, and their zero-inflated and zero-adjusted extensions models were considered for the count data, whereas the beta distribution and its inflated versions were evaluated for the proportions. Random effects were also included to consider the multilevel structure of the experiment. The zero adjusted negative binomial model has better fitted the count data, whereas the inflated beta distribution performed the best for the proportions. Both models allowed us to properly assess the effects of social separation, brushing, and genetic lineages on sheep behavioral. We may conclude that GAMLSS is a flexible framework to analyze animal behavior data.


Introduction
Animal behavior studies have supplied useful information on animal welfare in a wide range of situations that elicit different emotional states and are currently entwined in animal welfare understanding (BROOM and FRASER, 2015). Scientific sources of evidence of rich emotional capacities in animals have contributed to increasing concern and recognition that they are sentient beings (BROOM, 2014). Several approaches can be used to assess emotional states in animals. Recent studies have focused on behavioral responses, which are commonly used as inferences of emotional states in animals (BOISSY et al., 2007;BOISSY and ERHARD, 2014).
Such studies usually produce complex and unstructured data, recorded from audios, images and videos, among others. These data sources can yield a multiplicity of information, registered through a wide variety of random variables. For the sake of illustration, let us consider a specific animal behavior, such as a head movement or some eye reaction to some human stimulus. Based on records provided throughout the experiment, we may simply verify if each animal have presented or not some specific behavior, configuring a binary outcome; we may be interested in the number of times that such behavior was repeated, yielding a count variable; or yet in the proportion of time that it was expressed, given rise to a continuous rate; or in the amount of time before such behavior can be observed, a continuous variable subject to censoring, among several other possibilities.
In addition to the huge amount of variables of interest, animal behavior studies usually present several other factors that should be taken into account in the statistical analysis. For example, these studies are often designed such that the individuals (animals) are nested in groups, as litters or herds, yielding multilevel data (STEVENS et al., 2017;GRAHAM et al., 2018). Moreover, the case of longitudinal studies is common, where each animal is assessed in different time points (GÓRECKA-BRUZDA et al., 2017;VENTER et al., 2019), or spatial studies, where spatial coordinates are relevant for analyzing how the animals behave (PATTERSON et al., 2017;VILLEGAS-RIOS et al., 2017). It is well known that grouped, longitudinal and spatial data produce correlated observations, which should be properly analyzed in order to obtain valid inferences, as pointed out, for example, by Fitzmaurice et al. (2012). Other common constraints that may be verified in animal behavior data are, for example, nonlinear relationships; varying levels of dispersion, skewness and/or kurtosis; over (or under) dispersion; zero inflation (deflation); and missing or censored data. Therefore, the statistical analysis of animal behavior requires flexible models, able to deal with the aforementioned data features. For additional information on the constraints that frequently emerge in animal behavior studies we refer to Garamszegi (2016) and Valletta et al. (2017).
Generalized additive models for location, scale and shape (GAMLSS) configure a general framework for univariate regression models, known to be widely flexible due to the large number of available probability distributions, allowing to analyze data with different levels of skewness and kurtosis, zero inflation, mixed (continuous and discrete) behavior, among others. In GAMLSS we can model each distribution parameter by including covariates, random effects and smoothers, unlike other regression methodologies usually considered in animal behavior as, for example, linear and generalized linear models; linear and generalized linear mixed models; and generalized additive models, that only allow modeling a location (mean) parameter. By this way, several of the related constraints present in animal behavior data may be properly addressed. GAMLSS also allows the use of nonparametric as well as normal random effects models. Random effects are a suitable tool to deal with multilevel data and allow to properly accommodate the correlation structure resulting from repeated measure designs. Smoothers, on the other hand, are useful in modeling nonlinear relationships between the distribution parameters and continuous covariates. Some important smoothing additive term functions available in GAMLSS are based on P-splines, cubic splines, and local regression. For details regarding smoothing terms in GAMLSS, see Stasinopoulos et al. (2017).
There are several recent studies dealing with new developments and applications of GAMLSS. Debele et al. (2017) consider the implementation of GAMLSS to non-stationary flood frequency analysis, while De  describe Gaussian Markov random fields within a GAMLSS framework, and its application for spatial data analysis. Furthermore, the World Health Organization recommended GAMLSS for obtaining reference levels and centile curves for human populations in clinical research (see WHO, 2006 andWHO, 2007). Additional works on the GAMLSS methodology can be seen, for example, in Flatley et al.(2019); Nakamura et al. (2019);Smith et al. (2019). We strongly recommend the GAMLSS seminal book by Stasinopoulos et al. (2017), and the compendium of distributions currently available in GAMLSS (RIGBY et al., 2019).
This study aimed to explore the potential of GAMLSS methodology in animal behavior analysis. Therefore, based on data from a study on sheep behavior we seek to investigate whether social context, through the use of physical grids, influences behavioral responses of sheep submitted to brushing. We have also investigated whether social separation effects are influenced by emotional reactivity, by assessing genetically selected sheep for low versus high social motivation. Thus, 20 sheep belonging to two different genetic lineages were used. Amongst the main information resulting from this experiment we may point out two types of variables: the number of posture changes and the proportion (or percentage) of time that each posture was expressed. A variety of different postures were recorded, such as ear postures (raised-up, horizontal, and asymmetric); eyes posture (closed, half-closed and open); and feeding behavior (eating, ruminating, not eating or ruminating), among others.
For the sake of illustration, we have considered here one variable of each type: the number of ear posture changes and the proportion of time performing the horizontal ear posture. Amongst the behavioral indicators used to assess emotions in animals, ear behavior and ear posture changes have provided valuable data on how animal welfare can be improved (REEFMANN et al., 2009;BOISSY et al., 2011;TAMIOSO et al., 2017). Sheep submitted to positive events show high proportions of horizontal ear postures and fewer ear posture changes than sheep in negative situations. Therefore, through this research, it is expected to better understand sheep perception toward a presumed positive event, such as brushing, using ear postures as a behavioral indicator. A number of constraints we have pointed out in the analysis of behavioral data can be found in these data, such that they were used here as motivations of the usefulness of the GAMLSS framework in this context. It is worth noting that the proposed models can be easily extended to the other behavioral variables, such that analogous procedures can be adopted, for example, for the number of head posture changes, or the proportion of time that a particular body posture is expressed.
The remaining of this work is organized as follows. Section 2 describes the case study. Section 3 introduces the GAMLSS methodology. The analysis of count data, referring to the number of posture changes, is explained and the corresponding results are presented in Section 4, while Section 5 presents the analysis of the variable based on a proportion of time. Finally, Section 6 ends this work with our concluding remarks.

Case study
In this work we consider data from an experiment carried out at INRA experimental farm of La Fage, Roquefort, in the South of France (TAMIOSO et al., 2018). Twenty 15-month-old Romane female sheep of two different genetic lineages, selected according to their contrasted behavioral reactivity towards temporary social separation, low and high reactivity, were exposed to brushing by a familiar human (for more details, see Tamioso et al. (2020)). Ewes were also assessed pre-brushing, during brushing and post-brushing for a fixed period of time. In addition, testing was organized in three consecutive sessions: in sessions 1 and 3 one metal grid was used to separate the brushed animal from group members. In session 2, two identical metal grids 1.70 m apart was used to separate the brushed animal from group members. There was not any specific adaptation period for session 2.
In summary, there are three categorical covariates which are specified next. Lineage: a two-level factor that classifies animals as reactive (R+) or not reactive (R-) to temporary social isolation; Session: a three-level factor that defines the experimental session (session 1, session 2 or session 3); Moment : a three-level factor that defines whether or not the animal was under human intervention (pre, during or post brushing). The number of ear posture changes and the time performing the horizontal ear posture were registered in each moment and in each session. By this way, each animal contributed for the final data set with nine records, such that we have repeated measures from the same sheep.
The main goals of our study were to investigate the effect of temporary separation, by comparing the results across sessions; effect of phase related to the brushing procedure, by comparing the results in different moments across sessions; and effect of genetic lineages. Additionally, possible interaction effects were also of practical interest. The data sets used in this study, as well as the R scripts developed for the data analysis, are available online and can be assessed in http://www.leg.ufpr.br/doku.php/publications:papercompanions:gamlss.

Generalized additive models for location, scale and shape
To investigate animal behavior outcomes we have used generalized additive models for location, scale and shape (GAMLSS) implemented in the gamlss package (STASINOPOULOS et al., 2007) in R (R CORE TEAM (2017)). GAMLSS is a distributional regression approach that extends the well-known generalized linear models (GLMs) and generalized additive models (GAMs) such that all model parameters (e.g. location, scale and shape parameters) can be modeled as linear, nonlinear or smooth functions of covariates. A parametric (continuous, discrete or mixed) distribution is assumed for the response variable which may be heavy or light-tailed, and positively or negatively skewed.

GAMLSS framework
For the general definition of GAMLSS, consider independent observations Y i , i = 1, 2, . . . , n, that follow a probability (density) function f Y (y i |θ) parametrized by θ = (θ 1 , ..., θ p ) T , a vector of p up to four distribution parameters, denoted by (µ, σ, ν, τ ) T when p = 4. Usually, but not necessarily, µ and σ are the location and scale parameters while ν and τ refer to shape parameters. The GAMLSS is defined as where g k (·) is a known monotonic link function that relates the k-th distribution parameter to the predictor η k , β k = (β 1k , β 2k , . . . , β J ′ K k ) T is a parameter vector of dimension J ′ K , X k and Z jk are design matrices of order n × J ′ k and n × q jk , respectively. γ jk is a random variable of dimension q jk for which it is assumed γ jk ∼ N q jk (0, G −1 jk ), G −1 jk is the generalized inverse of the symmetric matrix G jk = G jk (λ jk ) of order q jk × q jk and λ jk is a vector of hyperparameters.

Parameter estimation
The regression parameters (β k ) and random effects (γ jk ) are usually estimated based on the maximization of the penalized log-likelihood (l p ), given by where l = n i=1 log(f (y i |θ)) is the log-likelihood function. The penalized loglikelihood reduces to the usual log-likelihood when there are not any random effects or smoothing terms in the model. The R package gamlss allows three options to maximize the log-likelihood presented in (2): (i) CG algorithm, (ii) RS algorithm and (iii) a combination of both methods denoted by mixed. The CG algorithm uses the first, second and cross-derivatives of the penalized log-likelihood with respect to the distribution parameters and performs better for distributions with potentially highly correlated parameters. The RS algorithm does not use cross-derivatives, is computationally less expensive, and is suitable for distributions with orthogonal parameters. In this paper we used the RS algorithm. For more information on the estimation process we refer to Stasinopoulos et al. (2017).

Specifying the probability distribution
For each response variable we have initially fitted different candidate probabilistic distributions. We compared the fitted models based on the Akaike Information Criterion (AIC) and considering the residual analysis, using the normalized (randomized) quantile residuals (DUNN and SMYTH, 1996), as suggested by Stasinopoulos et al. (2017). Other common measures available for comparing fitted models include the Generalized (Pseudo) R 2 and Generalized AIC (GAIC), which generalizes the AIC and BIC (Bayesian Information Criterion). See Stasinopoulos et al. (2017) for more details.
Quantile residuals are based on the idea of inverting the estimated distribution function for each observation to obtain random variables with uniform distribution, and next, normally distributed residuals. The main advantage of these type of residuals is that, whatever the distribution of the response variable, they always have a standard normal distribution when the assumed model is correct.
The normalized (randomized) quantile residuals are given byr is the cumulative distribution function of the standard normal distribution. Theû i 's are defined differently for continuous and discrete response variables. For continuous random variables we haveû i = F (y i |θ), that is, the fitted cumulative distribution function evaluated at y i , i = 1, 2, . . . , n. For discrete random variables, F (y|θ) is a step function with jumps at a set of integers y ∈ R Y and in this caseû i is defined as a random observation from a uniform distribution on [û i1 ,û i2 ] = F (y i − 1|θ), F (y i |θ) , for i = 1, 2, . . . , n.

Specifying the linear predictors
After we have selected a probabilistic distribution using the goodness-of-fit measures and graphical analysis of normalized (randomized) quantile residuals, we proceeded with the specification of the linear predictors. This step was based on the results of successive likelihood ratio tests (LRT), as we seek to compare nested models. Details on the specification of the linear predictors and the selection strategies are given in Sections 4 and 5.

Analysis of the ear posture changes
In this section we analyze the count ear posture changes. Posture changes configure an important animal behavior component, since a greater number of posture changes reflects a more agitated animal. Hence, higher counts are associated with less positive states (TAMIOSO et al. 2018(TAMIOSO et al. , 2020. According to Figure 1 we can observe that there is a remarkable strong right asymmetry noted in (a), which is an indication that there is overdispersion in the data and so the usual Poisson distribution may not be an adequate choice to model this data. Overdispersion, if not taken into account, can lead to underestimated standard errors and unreliable hypothesis tests (COX, 1983). The negative binomial distribution is usually considered to model overdispersion uncaptured by the Poisson distribution. We can see in (b) that the counts are usually lower for non-reactive animals and (c) shows that in session 1 there is a greater variability in response when compared to the other experimental sessions. Furthermore, we can observe in (d) that there is indication of differences in the frequency of changes in ear posture at different experimental moments. It may be noted, for example, that the counts are usually lower during brushing. We have also considered models able to accommodate a possible excess of zero counts in the response variable: the zero inflated Poisson, the zero inflated negative binomial, the zero adjusted Poisson and the zero adjusted negative binomial distributions. For zero inflated distributions, a discrete probability distribution, as the Poisson or negative binomial, has its probability at zero inflated, by adding a non null probability provided by a degenerate random variable at zero. Zero adjusted distributions, on the other hand, are mixed distributions that can be used when the non-zero part of the distribution is well fitted by some count model truncated at zero, whereas for the zero values, another probability distribution is assumed. Next we provide details on the probability distributions we have considered in this work.

• Poisson distribution
The probability function of the Poisson distribution, denoted by PO(µ), is given by where µ > 0. For the Poisson distribution we have E(Y ) = Var(Y ) = µ (equidispersion), which may be a very restrictive assumption when analyzing count data.
• Zero adjusted Poisson distribution where POtr(µ) is a Poisson distribution truncated at zero. Then Y has a zero adjusted Poisson distribution, denoted by ZAP(µ, ν), with probability function given by . The parameter µ is the mean of the Poisson component before truncation at zero and ν is the exact probability that Y = 0. The resulting marginal mean is E(Y ) = cµ while the marginal variance is Var(Y ) = cµ + cµ 2 − c 2 µ 2 .
• Zero adjusted negative binomial distribution Let Y = 0 with probability ν and Y = Y 0 with probability (1 − ν), where Y 0 ∼ NBItr(µ, σ) and NBItr(µ, σ) is a negative binomial truncated at zero. Then Y has a zero adjusted negative binomial distribution, denoted by ZANBI(µ, σ, ν), with probability function given by where µ > 0, σ > 0 and 0 < ν < 1, where Y 1 ∼ NBI(µ, σ). The parameters µ and σ represent the parameters of the negative binomial component before truncation at zero, respectively, while ν is the exact probability that Y = 0. By this way, the marginal mean is E For each of the aforementioned distributions, the regression model was specified with the fixed effects of session, moment and lineage, and all two-way interactions involving these three factors. Additionally, two random effects were included in the model: one at the animal level and another for animal within the sessions. These random effects are justified by the need to incorporate the non null correlations due to the repeated measures resulting from the multilevel design structure, such that each of the 20 animals contributes with 9 measures to the data set.
Let Y ijkl be the response variable for which i (i = 1, 2, . . . , 20) represents the animal, measured in session j (j = 1, 2, 3), at moment k (k = 1, 2, 3) and l corresponds to the lineage (l = 1, 2). The following general model specification was considered where µ ijkl is the location parameter of the conditional distribution of y, σ is a scale parameter and ν ikl is related to the excess of zeros. For all distributions, the initial regression structure for the location parameter was specified as where the parameters α (1) , β kl represent the corresponding two-way interaction effects. Further, u i is the animal random effect and v ik is the random effect for the animal nested within experimental session, for which we assume u i ∼ N(0, σ 2 U ) and v ik ∼ N(0, σ 2 V ), respectively. Additionally, for the zero inflated and zero adjusted distributions we have specified the following regression structure for the parameter corresponding to the zero portion We have additionally considered σ as a nuisance parameter. Special attention was devoted to the location and zero-excess parameters, due to their practical interpretations. It is possible that we could get a better fitted model if the covariates effects were also considered for σ. However, we have also took into account the risk of overfitting in our decision of considering σ constant. Table 1 shows the log-likelihood, degrees of freedom and AIC measures for the six fitted models. We may notice a poor performance for the Poisson distribution and its derived models that accommodate excess of zeros, which can be justified due to the presence of overdispersion in the data, which was not properly modeled by the Poisson distribution. On the other hand, the negative binomial distribution and its extensions presented better measures with the best performance observed for the ZINBI and ZANBI models. Although the ZINBI model has produced lower AIC than ZANBI, this model presented convergence problems and unreliable estimates and standard errors. For this reason, the ZANBI model was chosen. After we have selected the probability distribution, three nested models provided by two variations of the linear predictors were proposed and these were compared via LRT. The first model, denoted by Fit 1, is that specified by (4) and (5). In the second (Fit 2) we removed all two-way interactions from the linear predictor for µ, such that: and retained the predictor for ν as presented in (5). Finally, in the third model (Fit 3), besides the removal of the interaction effects in µ, as described in (6), we have also omitted all covariates in the linear predictor for ν: The results are given in Table 2, and we may verify that the second fit does not differ significantly from the original full model. However, a significant poorer fit is verified when the experimental effects were removed from the linear predictor of ν, in such a way that we chose the second fit as the final model.  Figure 2 shows the residual analysis for the chosen fit model. We may observe that the residuals present constant variance and no systematic pattern. In addition the QQ-plot and kernel plot of residuals do not exhibit any evident deviations from the standard normal distribution.   Table 3 presents the parameter estimates for the final model. Relative rates (RR) and odds ratios (OR) were obtained as exponentiated estimates for µ and ν parameters, respectively. Confidence intervals and Wald type hypothesis tests, based on the asymptotic normality of maximum likelihood estimators, were also provided. The results show that there was a significant effect of all variables under study considering a significance level of 5%, that is, the variations in ear posture change when the animals undergo intervention and social isolation, in addition to being different for animals of different lineages.
The results also indicate that the animals presented a lower rate of posture changes in sessions 2 and 3 than in session 1. It was still found a lower rate of changes in ear posture during or after brushing than before it. Note that the rate of ear posture changes is 0.35 times that observed before brushing, that is, we can expect about 3 times less changes in posture when the animals undergo intervention. It is also observed that the frequency of ear posture changes after the intervention is 0.65 times that observed before brushing. In addition, it was found that animals of reactive lineage moved 1.33 times more than non-reactive animals. The estimates referring to the inflation parameter can be interpreted as the propensity of the animal do not move its ears. It was found that animals have a higher chance of not moving their ears during brushing, a situation where animals are 15.8 times more likely to not move their ears when compared before brushing. A similar effect was found for experimental session for which there is a significant 14-fold increased chance of animals not moving their ears in session 3 as compared to session 1. In addition, the results indicate that animals belonging to the non-reactive lineage are more likely to remain with immobile ears than those considered reactive to temporary social isolation. The estimates for σ U and σ V are, respectively, 0.38 and 0.56. Finally, the LRT for comparison with a model including only fixed effects returned p < 0.001 reinforcing the need to take clustering into account in the analysis.

Analysis of the time expressing horizontal ear posture
Following, in this section we analyze the proportion of time that the animals remained with their ears in horizontal position. Horizontal ear posture is a sensitive indicator of positive states, such that a higher proportion of time expressing such posture evidences a more relaxed animal (TAMIOSO et al., 2020). Figure 3 shows in (a) that this response is considerably inflated at 0, such that in 72 occasions the animals did not remain with their ears in horizontal position at any time, which is equivalent to 40% of the collected measures, whereas there was only one case of proportion equals to 1. It is worth noting in (b) the greater variability in response for reactive animals and the same can be seen in (c) for measurements collected in session 1 when compared to the other experimental sessions. In addition, there are evidences of differences between the experimental moments (d). In order to model the proportion of time that the animals remained with their ears in horizontal position, the beta and inflated beta distributions were considered (see Ferrari and Cribari-Neto (2004), Ospina and Ferrari (2010) and Ospina and Ferrari (2012) for additional information on these models).

• Beta distribution
The beta distribution, denoted by BE(µ, σ), is a two-parameter distribution for a continuous random variable defined in the open interval (0, 1), and it is usually applied in modeling rates and proportions. The probability density function of the beta distribution is where α = µ(1 − σ 2 )/σ 2 and β = (1 − µ)(1 − σ 2 )/σ 2 , α > 0, β > 0 and hence 0 < µ < 1 and 0 < σ < 1. For the beta distribution we have E(Y ) = µ and Var(Y ) = σ 2 µ(1 − µ). Although the beta distribution is very flexible it has a drawback that it is not defined at 0 and 1. In order to use it in the presence of zeros and ones, we consider the following transformation suggested by Smithson and Verkuilen (2006): y c = [y(n − 1) + 0.5] /n where y is the observed response and n is the sample size.
• Inflated beta distribution The inflated beta distribution, denoted by BEINF(µ, σ, ν, τ ), is a four-parameter distribution for a continuous random variable defined in the interval [0, 1]. This distribution involves a mixture of three components: a mass probability p 0 at y = 0, a mass probability p 1 at y = 1 and a beta distribution defined in (0, 1) with probability (1 − p 0 − p 1 ). The probability density function of the inflated beta distribution is where W ∼ BE(µ, σ), p 0 = ν(1+ν+τ ), p 1 = τ /(1+ν+τ ) and p 2 = 1−p 0 −p 1 . Hence ν = p 0 /p 2 and τ = p 1 /p 2 . For the BEINF distribution we have E(Y ) = τ +µ 1+ν+τ and Var(Y ) = σ 2 µ(1−µ)+µ 2 +τ +(µ+τ ) 2 (1+ν+τ ) −1 1+ν+τ . The full regression model to analyze the proportion of time expressing horizontal ear posture for both the BE and BEINF distributions was specified as (1) kl represent the two-way interaction effects; u i is an animal specific random effect, and v ik is a random effect for the animal within the experimental session, for which we assume u i ∼ N(0, σ 2 U ) and v ik ∼ N(0, σ 2 V ), respectively. The regression structure for the parameter referring to the zero inflation for the BEINF distribution was specified as The parameter τ , referring to the one inflation, was not modeled through covariates effects since there was only one sample observation with Y = 1. The choice between the two distributions may be supported by the graphic analysis of the quantile residuals. The results presented in Figure 4 clearly show a poor fit of the beta distribution, for which the residuals exhibits a non normal, highly skewed distribution. On the other hand, residuals from the inflated beta distribution display satisfactory adherence to the normal distribution.
After we have chosen the inflated beta distribution we proceed by evaluating three nested regression models aiming to select the terms for the linear predictors. The first model, denoted by Fit 1, is defined by the regression structures in (7) and (8). In the second model (Fit 2) the interaction effects in the predictor for µ were removed, such that: and retained the predictor for ν as presented in (8). In the third model we have included no interaction terms for µ, as described in (9) nor covariates for ν, so that:    Based on the results presented in Table 4, we selected the original full model and proceeded with the interpretation of parameter estimates. Since we have used the logit link function in both linear predictors, we can interpret the exponentiated estimates as odds ratios. Once again, asymptotic CIs and hypothesis tests are given. Results are shown in Table 5. Additionally, Figure 5 shows the predicted mean values for each combination of factors for a better interpretation of the results.  We can notice that the proportion of the time the animals spent with the ears in a horizontal position are higher in session 1 than in the others, especially during brushing and in reactive animals. It was also noted that the lowest average proportions are observed in session 2, that is, in the one in which social isolation was imposed. The significant effect of the interaction between the second session and the moment indicates that in this section the brushing effect was less mild than in the other sessions, which can be seen in Figure 5. In addition, animals of reactive lineage remain longer with their ears in horizontal position than those belonging to the non-reactive lineage, mainly during human intervention and in the first experimental session.
Referring to the zero inflation, i.e., to the animals keeping their ears in some position other than horizontal, we can verify that it is higher under human intervention (during brushing), being approximately 3 times more likely than in the pre-brushing phase. We can also notice that in sessions 2 and 3 the animals are more likely to keep their ears in non horizontal postures than in session 1. The estimates obtained for σ U and σ V are, respectively, 0.0002 and 0.88, with p < 0.001.

Concluding remarks
In this paper we aimed to explore the flexibility of generalized additive models for location, scale and shape (GAMLSS) in modeling animal behavior data. Two types of random variables typically found in animal behavior research were analyzed, the first corresponding to a count variable (number of ear posture changes), and the other configuring a continuous proportion with inflation at zero and one (the proportion of time performing the horizontal ear posture). The main goals of this study were to investigate the effects of temporary separation, phase, and genetic lineage. The available data exhibited multilevel structure, overdispersion, zero and/or one inflation, and covariates effects in more than one model parameters. GAMLSS allowed us to properly analyze both data sets, considering all data properties and experimental settings, and providing reliable statistical results and consistent biological findings. The results of this study might help in devising strategies to promote positive welfare.
It is worth emphasizing that GAMLSS, the R package gamlss and its extensions offer several other possibilities for analyzing animal behavior data. They may be useful, for example, in fitting censored (interval) response variables (see the gamlss.cens package, by Stasinopoulos et al. (2018)), spatial data (gamlss.spatial package, by De Bastiani et al. (2018)), besides Bayesian (BAMLSS package, by Umlauf et al. (2018)) and boosting methods (gamboostLSS package, by Thomas et al. (2018)), among others. Moreover, it is well known that animal behavior variables are usually correlated, and these correlations cannot be taken into account through univariate regression models. Thus, multivariate GAMLSS and related methodologies, such as the multivariate covariance generalized linear models (BONAT, 2018), configure potential future research topics for animal behavior analysis. Alternative regression models for count and rate behavior variables could be also considered based on other distributions already implemented in the gamlss framework, as the double Poisson and Poisson inverse Gaussian (for counts) or simplex, logit normal or generalized beta type I (for continuous or mixed rates).