Ordinal data and residual analysis: Review and application

,


Introduction
In agricultural sciences, it is common to carry out experiments that result in polytomous data as a response of interest.These data assume values in a finite set of categories with nominal or ordinal scale (natural ordering between categories) and have a multinomial distribution regardless of this nature (Agresti, 2002).The models with the logit link function are the most used in the statistical analysis of these data.The proportional odds model (McCullagh, 1980) is widely used for the ordinal case with a smaller number of parameters due to the assumption of proportionality (Tutz, 2011).However, other alternatives can be considered, such as the cumulative probit model or the Proportional Hazards model with a complementary log-log link function (Agresti, 2010).When the proportionality assumption is not valid, the cumulative logit model (Williams & Grizzle, 1972) can be fitted to the data or the adjacent-categories logit model, for example (Ananth &Kleinbaum, 1997 andAgresti, 2002).Furthermore, one can assume another discrete multivariate distribution for the response variable, such as the Dirichlet distribution, which is the conjugate distribution of the multinomial in Bayesian inference (Ng et al., 2011).
When selecting a model, it is essential to assess the quality of its fit to the data as well as to validate its assumptions.The fitted model must describe the observed data well so as not to result in incorrect inferences.In this context, residual analysis plays an important role in detecting possible failures resulting from the fit and identifying outliers and/or influential points, becoming an integral part of any regression problem (Cook & Weisberg, 1982).McCullagh & Nelder (1989) paid substantial attention to defining residuals for Generalized Linear Models (GLMs), with Pearson and deviance residuals frequently used in the diagnostics of GLMs.However, these residuals do not apply to multinomial data due to the nature of the response variable.As the polytomous variable is multivariate, the ordinary residual given by the difference between the observed response and the estimated probability is a vector for each subject (Reiter & Kohnen, 2005).Therefore, diagnostic plots of residuals are difficult to interpret since their distribution is difficult to identify.Furthermore, few papers in the literature involve types of residuals that help validate models associated with polytomous data, and these are defined, in particular, for the case in which the response of subject results in only one of the categories.
For the ordinal case, Liu et al. (2009) presented the vector of cumulative residuals focusing on validating the proportional odds model with respect to the covariates of the linear predictor.However, it is not simple to interpret the behavior of these residuals in diagnostic plots, as is the case with residuals for continuous variables.Li & Shepherd (2012) and Liu & Zhang (2018) defined residuals that correspond to a single value per subject regardless of the number of categories.While the residual proposed by Li & Shepherd (2012) is obtained in the discrete space of the original response, in the approach used by Liu & Zhang (2018), a continuous variable replaces the ordinal variable, and the residual is defined through this new variable.Liu & Zhang (2018) compared the performance of the residuals so-called surrogate, with those proposed by Li & Shepherd (2012) in the residuals versus covariates plot and Quantile-Quantile plot (Q-Q plot) to assess the fit of the cumulative probit model with respect to mean structure, heteroscedasticity, and proportionality.The authors showed that the surrogate residuals presented expected behaviors in these plots for the model correctly specified to the data.In contrast, the residuals defined by Li & Shepherd (2012) showed unusual patterns that did not allow concluding in favor of the correct model.
The aim of this work is to present a review of models and residuals for polytomous ordinal data, considering the relevance and need for studies and research in this area.As a specific case, we show the performance of the surrogate residuals to evaluate the cumulative logit model for ordinal response.As a motivational study and application, it is presented the research carried out with Tambaqui fish (Colossoma macropomum), in which a type of histopathological alteration was observed in the liver fish.Therefore, in this study, the response variable is the severity of lesion found in the fish liver (natural ordering), which was classified as mild, moderate, and irreversible.Also, it is verified the relationship of the classifications with the different gene expressions of the Tambaquis.This species is a source of aquatic protein widely consumed in the North region of Brazil and has attracted significant interest from fish farmers from other countries (Lopes et al., 2016).Given the large production of Tambaqui in the country, the aquatic environment and the management of these fish must be appropriately controlled to generate a healthy population, not causing losses in productivity (Correa et al., 2018).

Models for ordinal response
The multinomial distribution is the most important and usual one for random variables associated with categorical data.Thus, it is the distribution that is assumed, except for overdispersion, in classic models with polytomous responses (Agresti, 2002).Let it be a multinomial trial, that is, an experiment that admits J possible and mutually exclusive outcomes, whose probabilities are denoted by π 1 , π 2 , ..., π J such that 0 ≤ π j ≤ 1, j = 1, 2, ..., J , and Consider m identical and independent trials, which means that the probabilities of occurrence of the results are constant for each trial and that the result obtained in one trial does not interfere with the result of the other.The components of the random vector Y = (Y 1 , . . ., Y J ) ′ give the cell counts in categories 1, . . ., J .Then, the random vector follows a multinomial distribution with parameters m and π = (π 1 , . . ., π J ) ′ , Y ∼ Multi(m, π), and probability mass function given by where y j ∈ {0, 1, . . ., m} and For the category j the result y j has mean and variance given by E(Y j ) = mπ j e Var(Y j ) = mπ j (1 -π j ), respectively.Furthermore, the covariance between y j and y j ′ , ∀j ̸ = j ′ , j ′ = 1, . . ., J , is obtained by Cov(Y j , Y j ′ ) = -mπ j π j ′ , and that the marginal distribution of each y j is binomial.

Cumulative logit model
The cumulative logit model (Williams & Grizzle, 1972) is a multivariate extension in the class of generalized linear models used to model the dependence of an ordinal response on discrete or continuous covariates.In this context, the response variable Y i takes on a value in the set 1, 2, . . ., J for the i-th subject, i = 1, 2, . . ., n, with the ordered categories 1 < 2 < . . .< J and following the multinomial distribution.Then, the cumulative logit models with the canonical link function can be used to describe the functional relationship between the response and covariates of the study.According to Agresti (2010), models that consider the natural order of the response can produce more powerful results than models that ignore ordinality.This model is defined by: where x i = (x i1 , x i2 , . . ., x ip ) ′ is the vector of the p covariates for the i-th subject, β j = (β j1 , β j2 , . . ., β jp ) ′ represents the vector of regression parameters and α j is the intercept, with j = 1, 2, . . ., J -1, these ones constitute the systematic part of model.Here for the multinomial response Y i , γ ij (x i ) is the cumulative probability of the subject i until the j-th category, that is, ) the probability of the (marginal) response in the j-th category, more precisely, In the cumulative logit model, the regression parameters are not constant for the j logits, i.e., β j can vary according to each response category.The estimation of the parameters of the model ( 1) is generally performed using the maximum likelihood method, whose likelihood function for the random sample of size n is given by , where y ij = 1 if the response of subject i, i = 1, . . ., n, belongs to the category j, j = 1, . . ., J , is the vector with the parameters to be estimated.It is necessary to use iterative methods such as the Newton-Raphson method to maximize L and obtain the maximum likelihood estimators of the parameters (Agresti, 2002).Additionally, as is a classic in generalized linear models, under conditions of regularity, in the ordinal case, asymptotically θ has a normal distribution, that is, θ ∼ N(θ, ℑ -1 ), where ℑ is the Fisher information matrix.This matrix is of fundamental importance in the construction of hypothesis tests and asymptotic confidence intervals for the elements of θ (via the Wald method), since the square root of the main diagonal elements are the standard errors of the estimators.More details can be found at McCullagh (1980).
An alternative to model (1) is the proportional odds model, which assumes that the effects of the covariates are the same for each logit j, resulting in a more parsimonious model, that is, with a smaller number of parameters (Bilder & Loughin, 2014).The proportional odds assumption results in the simplest fit with easy interpretation, but it should always be carefully verified (Lemos et al., 2015).

Proportional odds model
The simplest model in the class of cumulative logit models involves parallel regressions on the ordinal scale and assumes equivalent proportions by assuming the same regression parameter for all categories.This model, called the proportional odds model, was introduced by McCullagh (1980) and can be expressed by where x i = (x i1 , x i2 , . . ., x ip ) ′ is the vector of covariates for the subject i, β = (β 1 , β 2 , . . ., β p ) ′ represents the vector of regression parameters, α j is the intercept, and the last category J as the reference.Here for the multinomial response Y i , γ ij (x i ) is the cumulative probability of subject i until the j-th category, that is, . ., J .The probabilities π ij (x i ) are obtained for the model ( 2) by means of subtractions given by where As the effects of the covariates are equal, the model assumes that the effects on the logit are identical for all categories of the response variable.Then, the J -1 logits are shifted only as a function of the intercept (Bilder & Loughin, 2014).According to Agresti (2007), the maximum likelihood procedure can be used to estimate the parameters of the model (2), with a likelihood function for the random sample of dimension n described by , where y ij = 1 if the response of subject i, i = 1, . . ., n, belongs to category j e y ij = 0 otherwise, j = 1, . . ., J , with J j=1 y ij = 1 and θ = α 1 , . . ., α J -1 , β ′ representing the vector of parameters.
According to McCullagh (1980), the Newton-Raphson method with Fisher scoring can be used to obtain parameter estimates, converging rapidly even with poor initial values.
The odds ratio is a very intuitive and used way to interpret the parameters estimated by the proportional odds model.Consider two subpopulations characterized by vectors x 1 and x 2 , then the cumulative odds ratio for the two subpopulations is given by where the odds of occurring As stated in Bilder & Loughin (2014), the cumulative odds ratio remains the same regardless of the category j used, and this is due to the assumption that the effects of the covariates are the same for all categories.
As the proportional odds model is a particular case of the model (1), the proportionality assumption can be verified through the likelihood ratio test (LRT) with the following hypotheses and with the statistic of the test given by where L H 0 is the likelihood function under the null hypothesis H 0 , i.e., referring to model (2) and L H 1 is the likelihood function under the alternative hypothesis H 1 , i.e., referring to model (1).Here, Λ follows an approximate Chi-square distribution, in which the degrees of freedom, m, are obtained by the difference between the numbers of the parameters under the hypotheses H 0 and H 1 .If the null hypothesis is not rejected at the 5% significance level, then the proportional odds model can be fitted to the data (Lemos et al., 2015 andGiolo, 2017).
The proportionality assumption can be verified in two ways: global and subject.Globally, all model covariates are considered, while subjectly, it is considered covariate by covariate.In the case of rejection of the null hypothesis for part of the covariates, that is, some covariates have the property of proportional odds and others do not, an alternative is the partial proportional odds model (Agresti, 2010).

Partial proportional odds model
The proportional odds assumption is not always achieved in practice.A model proposed by Peterson & Harrell Jr (1990), an extension of the proportional odds model, can be used when part of the covariates violates this assumption.
Consider the vector x i with the values of p covariates for the i-th subject that present proportional odds and the vector z i with the values of q (q ≤ p) covariates that do not, so the partial proportional odds model is given by where β = (β 1 , β 2 , . . ., β p ) ′ and ρ j = (ρ j1 , ρ j2 , . . ., ρ jq ) ′ sare the vectors of regression parameters, α j is the intercept and the last category taken as a reference.Here, the vector ρ j describes the effect of non-proportionality for each j-th cumulative logit associated with the vector z i .In this model, J -1 intercepts, p coefficients referring to the vector β, which are independent of the compared categories, and q(J -1) coefficients referring to the vector ρ j are estimated.Furthermore, . ., J , and the probabilities π ij (x i , z i ) for th model ( 3) are obtained in an analogous way to those obtained for models (1) and ( 2), so , where The estimation of parameters can also be performed using the maximum likelihood method for the random sample of size n (Agresti, 2010).Considering y ij = 1 if the response of subject i, i = 1, . . ., n, belongs the category j, j = 1, . . ., J , y ij = 0 otherwise and J i=1 y ij = 1, the estimators of the model (3) can be obtained by maximizing the likelihood function (or its logarithm) given by corresponds to the vector of parameters to be estimated.
The estimates can be obtained using the step-halving technique in the modified Gauss-Newton algorithm that ensure, in each iteration, an increase in the likelihood logarithm (Peterson & Harrell Jr, 1990).
The adjacent-categories logit model is also an alternative when the proportionality assumption is not satisfied.It considers the ratio between the probabilities of successive categories rather than the cumulative probabilities.We reinforce that all models mentioned here are based on the multinomial distribution for the response variable, i.e., a random part of the model.They differ in the structure of the linear predictor, the systematic part of the model.Additionally, it is possible to find this model and others for ordinal data in Ananth & Kleinbaum (1997), Agresti (2002), Agresti (2007), Agresti (2010), Tutz (2011), Bilder & Loughin (2014), Giolo (2017), among others.

Residuals for ordinal data
After fitting a model to the data, it is essential to verify whether its assumptions are satisfied and identify subjects that may disproportionately interfere with the results obtained.Through an analysis of the residuals, it is possible to study the robustness of the fitted model in terms of the various aspects that involve its formulation and the estimates of its parameters, detecting potential problems, and improving the fitting process (Souza, 2006).

Ordinary Residual
For the class of models with a polytomous categorical response, the ordinary residual associated with the i-th subject, i = 1, . . ., n, is a vector J × 1 defined by (Reiter & Kohnen, 2005) where y i = (y i1 , y i2 , . . ., y iJ ) ′ is the observed vector with y ij = 1 if the response of the subject i belongs to the category j and y ij = 0 otherwise, πi = ( πi1 , πi2 , . . ., πiJ ) ′ is the estimated probabilities vector.The only positive element in this vector pertains to the observed outcome for the subject.This vector may not be informative in the diagnostic techniques for analyzing residuals since its asymptotic distribution is unknown.

Cumulative Residual
Specifically for the proportional odds model, defined in section 2.2, Liu et al. (2009) presented the cumulative residuals for a binary response (by collapsing the categories) and the vector of cumulative residuals considering the original response.For the multivariate case, the vector of cumulative residuals, J × 1, for each subject is expressed by where ′ is the vector of cumulative probabilities for the i-th subject.The authors used the sum of this residual vector in graphical and numerical methods to assess the goodness-of-fit of the model.The methods generalize those developed by Arbogast & Lin (2005) for the logistic regression model with binary responses.However, diagnostic plots associated with residuals are difficult to interpret.

LS Residual
Considering the models that assume the assumption of proportionality for the regression parameters, Li & Shepherd (2012) proposed a residual that is a single value per subject, regardless of the number of ordered categories.This residual, called LS, is obtained by the difference between two cumulative probabilities, and the authors examined several properties to apply it to the available diagnostic tools.The residual associated with a subject considering the model 2 is obtained by with its value varying in the numeric interval of [-1, 1].The Q-Q plot of this residual is obtained compared to the theoretical quantiles of a Uniform distribution in [-1, 1].However, the residual is defined on the discrete space of the response variable, and its conditional distribution can vary according to the covariates.These facts make it difficult to analyze the residuals in diagnostic plots since they do not produce the expected patterns.According to Liu & Zhang (2018), the use of this residual is limited to verifying its zero mean under the correct model.

Surrogate Residual
The residual defined by Liu & Zhang (2018) is also a single value per subject for the models that assumes the proportional odds.Consider the model (2) and a latent variable given by Z i = -β'x i + ε i , i = 1, 2, . . ., n, where ε 1 , . . ., ε n is a random sample of the variable ε which follows a standard logistic distribution, ε ∼ log(0, 1), with probability density function and cumulative distribution function, respectively, given by where u ∈ R. The mean and variance of ε are E (ε) = 0 and Var (ε) = π 2 3 , respectively.The concept of latent variable induces a joint distribution of the variables Y i and Z i determined by Thus, the marginal distribution of the ordinal variable Y i is the same as the distribution specified by the model (2).The authors defined a continuous variable S i based on the conditional distribution of Z i given Y i , i.e., S i follows a truncated distribution of Z i in the interval α j-1 ; α j given Y i = j.Therefore, the surrogate residual is obtained by the difference between the surrogate variable and its expected value, with the expression given by where E 0 (.) and E(.) denote, respectively, the mean of variables S i and Z i .If the model ( 2) is specified correctly, the variable S i follows the same distribution of Z i and the residual R S i , which is also a continuous variable, has the following properties: 3 , a constant does not depend on x i ; iii) Reference distribution: Independent of x i , the empirical distribution of R S i approximates of the standard logistic distribution, that is, R S i ∼ G (.).These properties allow an analysis of residuals in practically all existing diagnostic tools for continuous variables Liu & Zhang (2018).As the residuals are obtained by random sampling, diagnostic plots may vary from one sample to another (especially for small samples).The authors presented a bootstrap algorithm for the residual (5) similar to the bootstrap algorithm used in linear regression proposed by Efron (1979)  2) Using the bootstrap sample obtained in step 1, perform the conditional sampling procedure presented in this section to generate a sample of the surrogate residuals given by R S * 1b , R S * 2b , . . .R S * nb .Thus, it is possible to examine the discrepancy between the empirical bootstrap distributions and the reference distribution (standard logistic).As the bootstrap samples are drawn independently, the behavior of B × n surrogate residuals is examined in the plot of residuals versus covariate (or fitted values), while the median of the B bootstrap distributions is examined in the Q-Q plot.

Diagnostic techniques
Several diagnostic techniques based on residual analysis can assess the goodness-of-fit of a statistical model.These can be informal through residual plots or formal when using tests.The tests provide a p-value referring to a tested hypothesis.At the same time, the graphical representation is an important exploratory diagnostic feature that can reveal which components of the model were not correctly specified.
When fitting a linear regression model, the Shapiro-Wilk test (Shapiro & Wilk, 1965) is generally used to verify the normality assumption of residuals.On the other hand, the Kolmogorov-Smirnov test (Kolmogorov, 1933) is a widely known test that considers continuous models other than the linear regression model.Through this test, it is possible to examine the degree of agreement between the empirical distribution function of the residuals concerning the theoretical distribution function of reference (Dufour et al., 1998).In addition, a simple way to visualize the shape of the residual distribution is through a histogram, making it possible to compare the result obtained with the shape of the normal distribution or any other distribution.
Consider R 1 , R 2 , . . ., R n a random sample of residuals with empirical distribution function and G(c) the theoretical distribution function of reference.The hypotheses of the Kolmogorov-Smirnov test are given by and test statistic between the two distribution functions.For a significance level α = 5%, the H 0 is rejected if the statistic T n exceeds the quantile value of 1α as given by the table of quantiles for the Kolmogorov test statistic.In case of non-rejection of the null hypothesis, R 1 , R 2 , . . ., R n is a random sample from the theoretical distribution function.
Although goodness-of-fit tests provide a p-value that indicates how strong the evidence (observed data) is against the null hypothesis, they may fail in certain circumstances, for example, when the sample size is small.Generally, graphical techniques can be more informative, providing a better diagnostics of model adequacy than hypothesis testing (Moral et al., 2017).Among the different types of diagnostic plots, some principals are (Paula, 2013 i) Residuals versus covariates: indicates whether the systematic part was incorrectly specified, with the need to include higher-order terms or transform the quantitative covariates into the linear predictor.The expected pattern of this plot is a zero-centered distribution of residuals with constant amplitude; ii) Residuals versus fitted values: the behavior of the residuals in this plot must be the same as described in item (i) for a well-fitted model.This plot can reveal the existence of heterogeneity of variance in addition to outliers; iii) Normal and half-normal plots: they are widely used for the diagnostics of the model, being possible to detect outliers and identify failures in the specification of the link function or distribution of the random component.The residuals should follow approximately a straight line with a slope of 45º for a well-fitted model.
Under the normality assumption, the normal plot of the residuals against the expected sorted values of the standard normal distribution, which is approximated by while in the half-normal plot, the absolute values of the residuals (even with unknown distribution) are compared concerning the expected order statistics of the half-normal distribution, obtained by where Φ -1 is the standard normal distribution function, with i = 1, . . ., n and n corresponding to the sample size.However, the behavior interpretation of the points in these plots can be subjective, and it is difficult to point out other causes for unavoidable irregularities.To assist in visual analysis, Atkinson (1985) proposed adding a simulated envelope to these plots.So, it is possible to observe the proportion of points randomly distributed within the envelope and decide whether the observed residuals are consistent with the fitted model.The envelope is simulated with a confidence level of 95% that contains the residuals, i.e., there is evidence of a good fit when the number of points outside the envelope is below or equal to 5%.
In addition, a sensitivity analysis based on a set of tools (such as leverage, case-deletion, or local influence analysis) can be designed to evaluate changes in the fitted model when some perturbation is imposed on the data or assumptions of the model (Singer et al., 2017).This will be briefly presented below given that the focus of this paper is on the residual analysis.
It is important to examine the existence of one or more points poorly fitted by the model (do not follow the same pattern as the others) and may cause a significant impact on some characteristics of interest, such as the parameter estimate or the corresponding standard error (Singer et al., 2017).A simple technique introduced by Cook (Cook, 1977) that can be used is the deletion, which measures the impact on the fit of the model by considering all the subjects with the fit when deleting a particular subject from the sample.Consider θ and θ(i) the estimated maximum likelihood vectors from the sample with all subjects and the sample without subject i, respectively.An indicator of the influence of i-th subject can be calculated by θθ(i) .If the estimates differ substantially, the subject can be considered influential.
A measure that also can be used to assess the influence of the i-th subject is the likelihood displacement (Cook & Weisberg, 1982).This measure verifies the distance between the two likelihoods, being given by where l( θ) and l( θ(i) ) are, respectively, the likelihood logarithms of the parameters obtained from the sample with all points and the sample without the i-th subject.When it is not possible to obtain an analytical form for LD i , a quadratic approximation by Taylor series leads to the following result: where is the Fisher information matrix, which is estimated by substituting θ.Generally, it is not possible to obtain a closed form for θ(i) , and the one-step approximation is used that takes the first iteration of the iterative process by the Fisher score method when it starts at θ (Paula, 2013).
Another indication of an influential observation is through leverage point.The leverage value measures the potential for a subject to have a large effect on the fitted regression line, being defined as a measure of how far a particular case is (based on only predictor values) from the average of all cases.Also, looking at residuals doesn't help to detect leverage points since they don't necessarily fall off the regression surface (Simonoff, 2003).
Finally, there are several techniques that can help in the diagnostics of the model, which are described in several papers such as Junior & Veiga (2020) that assessed the local influence and the likelihood displacement measure for diagnostics in normal and logistic regression models.Details about these diagnostic measures are covered in Cook & Weisberg (1982), McCullagh & Nelder (1989), Turkman & Silva (2000), among others.It is highlighted that a point should only be excluded as a last alternative after several attempts to accommodate it in the fit, such as through transformations or including covariates (Silva, 2003).

Material
As an application, it is considered the data from the experiment conducted by Marques (2018) regarding the histopathological alterations found in the livers of Tambaqui fish (Colossoma macropomum) at the Biofish-Aquicultura farm based in Porto Velho-RO from January 2015 to October 2016.In this experimental study, juvenile fish were anesthetized and marked using a microchip in the ventral portion (Figure 1), applying Methylene Blue to the inserted site to prevent infection.After recovering from anesthesia, the Tambaquis were managed in an excavated pond with approximately 600m 2 of water, where they received the same food three times a day.In the end, the fish were fasted for 24 hours, collected with a trawl, and anesthetized when transported to water tanks for slaughter.The pituitary gland was collected for gene expression analysis, placed in a stabilizing solution (RNAlater), and stored at -80ºC until the moment of RNA extraction.With the DNA Analyzer 4300 equipment, two different types of genotypes, 122 and 130, were obtained.Small organ fragments were collected and properly stored for the liver histopathology analysis at the Laboratory of Ecology of Reproduction and Recruitment of Marine Organisms, Oceanographic Institut, USP/SP.The histopathological alterations were photomicrographed, Figure 2, and ordered according to the severity of the lesions, being classified as mild, moderate, and irreversible.Images of the lesions were obtained using the AXIOSKOP-ZEIS photomicroscope.The author made available 21 data from fish with genotype 122 and 21 from fish with genotype 130, totaling a sample of size equal to 42, in which was verified the relationship between the severity of lesions with the different gene expressions of Tambaqui.According to Marques (2018), the liver needs to function properly for a healthy fish population.

Methods
In this application, the response represents the histopathological alteration obtained in the liver of the fish associated with a different genotype, and the degree of severity of the lesions (from less to more severe) depends on this classification.Then, the response variable Y i , i = 1, 2, . . ., 42, has an ordinal scale assuming values in the set {1, 2, 3}, i.e., Y i = j represents the response of the i-subject in the category j, j = 1, 2, 3, where 1-mild, 2-moderate, 3-irreversible with 1 < 2 < 3.In this context, the corresponding observed vector is y i = (y i1 , y i2 , y i3 ) ′ , where y ij = 1 if the response referring to the fish i belongs to the category j and y ij = 0, otherwise.The genotype covariate is a factor, being incorporated into the model through the dummy variable.First, to test the proportionality, the likelihood ratio test described in section 2.2 will be used considering the model with the main effect given by where α j is the intercept, β j is the parameter associated with the genotype effect on the j-th logit.
Here, the third category is used as a reference.Using the standard parameterization, x i = 0 for the i-th fish with genotype 122 and x i = 1 for fish i with genotype 130.
If the proportionality condition is not violated, proportional odds are assumed.Otherwise, the model ( 6) is used to proceed with selecting the linear predictor.Under the proportionality assumption, the sequential proportional odds models are expressed by Model 1 -Null: Model 2 -Genotype effect: The likelihood ratio test (LRT) is used to select the structure of the linear predictor, verifying if there is an effect of genotype in the classification of severity found in the Tambaqui liver, that is, if H 0 : β = 0 is true or false.The test statistic is given by where l H 0 ( α) is the logarithm of the null model likelihood function and l H 1 ( α, β) is the logarithm of likelihood function of the model with genotype effect, with expressions given by and with α = (α 1 , α2 ) ′ .The estimates of the parameters of models ( 1) and ( 2) are obtained by the maximum likelihood procedure as described in the review chapter, section 2.2.The null model has only the intercept effect (2 parameters), and model 2 takes into account the intercept and genotype effect (3 parameters) under the null hypothesis has Λ ∼ χ 2 1 .Once the genotype effect is significant, confidence intervals (CIs) are constructed for the estimated probabilities for each response category and comparisons between observed and estimated proportions.In this way, simultaneous confidence intervals of 100(1α)% are given by (see May & Johnson, 1997) where χ 2 (α,l) is the point from a chi-square distribution with l = J -1 = 2 degrees of freedom and α = 0, 05 is the significance level.The estimated probabilities are expressed by Next step, for the fitted model validation, the surrogate residuals are used as described in section 3.4.Thus, with the data and the model, the conditional distribution of Z i ∈ (α j-1 ; αj ) given Y i = j is obtained by substituting the parameter estimates αj 's and β where the latent variable is Z i = -βx i +ε i and ε i ∼ Log(0, 1).A random sample s i , i = 1, 2, . . ., 42, is obtained from this distribution, and the i-th surrogate residual is given by Once obtained the residuals, it is possible to compare their empirical distribution function graphically with the standard logistic distribution function.Also, the bootstrap algorithm described in section 3.4 is used with 10 replications because of the sample size.The informal and formal techniques to evaluate the residual performance are the following: a) histogram, b) half-normal plot, c) the plot of residuals versus covariates, and c) the Kolmogorov-Smirnov test as described in section 4.
The analysis and estimation of model parameters were performed by the clm(.)function of the ordinal package (Christensen, 2013) and the resids(.)function of the sure package (Greenwell et al., 2018) to obtain the surrogate residuals.The ks.test(.)function of the dgof package (Arnold & Emerson, 2011) was used to obtain the p-value of the Kolmogorov Smirnov test.Finally, the hnp(.)function is used for the half-normal plot with a simulated envelope, implemented in the hnp package (Moral et al., 2017).All are available in the R software (R Core Team, 2020).

Results and Discussion
Initially, an exploratory analysis was carried out to describe the fish data set.The frequencies of mild, moderate, and irreversible lesions were obtained for each type of genotype (Figure 3), in which one can observe the differences according to classifications.The liver alteration classified as irreversible had a higher frequency in fish with genotype 122 than in fish with genotype 130. On the other hand, fish with genotype 130 had higher frequencies of mild and moderate lesions than fish with genotype 122.Then the cumulative logit and proportional odds models were fitted to test proportionality.It was verified evidence in favor of the proportional odds model by the LRT (pvalue = 0.8667).Afterward, the sequential proportional odds models were fitted and compared using the LRT as well.The model that considers the genotype effect was selected (p-value= 0.02714).
Based on this result, it is concluded that the type of genotype contributes to explaining the lesion classification in the liver of the Tambaqui fish in the study carried out by Marques (2018).
The estimated parameters and standard errors for the model with genotype effect are presented in Table 1.The expressions in terms of the cumulative logits for the proportional odds model with genotype effect are expressed by log γ 1 (x) 1γ 1 (x) = -3.1289+ 1.3779 x and log γ 2 (x) 1γ 2 (x) = -0.9079+ 1.3779 x.
The interpretation of the estimated parameter is generally performed through the odds ratios.The estimate of the genotype effect parameter is 1.3779 (Table 1), which indicates a tendency towards classification in the less severe categories in fish with genotype 130, as observed in the exploratory analysis.Therefore, the odds of the lesion being classified as mild (in relation to moderate or irreversible) in fish with genotype 130 was approximately 3.97 times the odds of being classified in fish with genotype 122.The same conclusions can be obtained considering the odds of the lesion being classified as mild or moderate in relation to irreversible, which occurs due to the proportionality assumption assumed by the model.
The predicted probabilities for each response category in the different types of genotype, with their respective confidence intervals, are presented in Table 2. Fish with genotype 122 showed irreversible liver alteration with a probability of 71.26%, while for fish with genotype 130, this occurs with a probability of 38.46%.Therefore, fish with genotype 122 tend to have more severe liver lesions than fish with genotype 130.As shown in Table 2, the confidence interval has greater amplitude due to the relatively small sample size.The observed and estimated proportions by genotype can be seen in Figure 4. Visually, the values are close to each other, showing that the proportional odds model includes the genotype effect is reasonable for describing the proportions of lesions observed in the study conducted by Marques (2018).The validation of the model assumptions was verified by the surrogate residuals analysis using bootstrap replications due to the sample size.Observing the histogram, Figure 5, the residual distribution presented a shape similar to the standard logistic distribution (red line), which is symmetrical, similar to the normal distribution but with heavier tails.The values for mean and variance were approximately 0.002 and 3.176, respectively.Furthermore, the p-value of the Kolmogorov-Smirnov test was approximately 0.729, which indicates in favor of the hypothesis that the surrogate residuals follow a standard logistic distribution.The half-normal plot with a simulated envelope for the surrogate residuals was presented in Figure 6.There is evidence that the observed data are a plausible realization of the fitted model since no systematic deviation pattern was observed with all the points inside the envelope.Thus, the model with the genotype effect can be used to analyze the data.Half-normal plot with a simulated envelope (the default is 0.95 confidence level) for the surrogate residuals to assess the fit of the model with genotype effect in the study of Marques (2018).
As in this model, a covariate is a factor, using the plot of residuals versus covariate is inappropriate.The boxplot of residuals was obtained for each genotype (Figure 7), which revealed medians of residuals close to zero.In addition, the residual distributions present symmetrical tendency, similar variability, and the presence of outliers per genotype.The large residuals in the Figure 7 refer to subject #32 for genotype 122 and to subjects #8, #13, and #15 for genotype 130.The model was fitted without these subjects to assess the impact on the estimates of model parameters.The parameter estimates and the related standard errors, in parentheses, are shown in Table 3.The variations between the estimated parameters (and the standard errors) were not disproportionate with the exclusion of subjects by genotype from the sample (Table 3), indicating that these points do not have a high influence on the fit.Thus, the entire inference based on the complete sample remains valid, and the choice of another model could lead to inadequate conclusions.Finally, the results were satisfactory, contributing to the validation of the model that provided a good fit for the data.

Figure 1 .
Figure 1.Microchip inserted in the juvenile of Tambaqui in a study carried out by Marques (2018) at the Biofish-Aquicultura farm.

Figure 2 .
Figure 2. Morphology of the liver tissue of the Tambaqui fish with the histopathological alterations in the experiment carried out by Marques (2018) at the Biofish-Aquicultura farm.A-Ductal hypertrophy (black arrows); B-Hemosiderosis; C-Cholestasis (black arrows); D-Focal necrosis (blue arrow) and Congestion of vessels and sinusoids (black arrows).

Figure 3 .
Figure 3. Frequencies of mild, moderate, and irreversible lesions in the liver of Tambaquis by type of genotype (122 and 130) in the study carried out by Marques (2018) at the Biofish-Aquicultura farm.

Figure 4 .
Figure 4. Observed proportions for the mild, moderate, and irreversible lesions and proportions estimated by proportional odds model with genotype effect in the study of Marques (2018).

Figure 5 .
Figure 5.Histogram of surrogate residuals related to the proportional odds fitted model (genotype effect) to the fish data in the study ofMarques (2018)

Figure 6 .
Figure 6.Half-normal plot with a simulated envelope (the default is 0.95 confidence level) for the surrogate residuals to assess the fit of the model with genotype effect in the study ofMarques (2018).

Figure 7 .
Figure 7. Boxplot of surrogate residuals per genotype to assess the proportional odds fitted model to the fish data in the study of Marques (2018).
to account for the variability of conditional sampling.It consists of repeatedly resampling the observed data, generating new data sets, and finding characteristics of interest in the population studied.The algorithm for obtaining the b-th bootstrap replication of surrogate residuals, b = 1, 2, . . ., B, is given in two steps (Liu & Zhang, 2018): 1) Generate a bootstrap sample of size n through sampling with replacement of the original data * 2b , . . ., x * nb , Y * nb .

Table 1 .
Marques (2018)ession parameters of the proportional odds model with the effect of genotype selected for analysis Tambaqui in a study carried out byMarques (2018)

Table 2 .
Estimated probabilities and 95% confidence intervals (in parentheses) in each of the response categories for fish with genotypes 122 and 130 were obtained by fitting the proportional odds model with the genotype effect

Table 3 .
Estimated Parameters of proportional odds model with genotype effect by excluding the subject #32 for genotype 122 and the subjects #8, #13 and #15 for genotype 130 the fish data