Bayesian Linear Regression and Natural Logarithmic Correction for Digital Image-Based Extraction of Linear and Tridimensional Zoometrics in Dromedary Camels

This study evaluates a method to accurately, repeatably, and reliably extract camel zoometric data (linear and tridimensional) from 2D digital images. Thirty zoometric measures, including linear and tridimensional (perimeters and girths) variables, were collected on-ﬁeld with a non-elastic measuring tape. A scaled reference was used to extract measurement from images. For girths and perimeters, semimajor and semiminor axes were mathematically estimated with the function of the perimeter of an ellipse. On-ﬁeld measurements’ direct translation was determined when Cronbach’s alpha (C α ) > 0.600 was met (ﬁrst round). If not, Bayesian regression corrections were applied using live body weight and the particular digital zoometric measurement as regressors (except for foot perimeter) (second round). Last, if a certain zoometric trait still did not meet such a criterion, its natural logarithm was added (third round). Acceptable method translation consistency was reached for all the measurements after three correction rounds (C α = 0.654 to 0.997, p < 0.0001). Afterwards, Bayesian regression corrected equations were issued. This research helps to evaluate individual conformation in a reliable contactless manner through the extraction of linear and tridimensional measures from images in dromedary camels. This is the ﬁrst study to develop and correct the routinely ignored evaluation of tridimensional zoometrics from digital images in animals.


Introduction
Zoometry, or the measurement and comparison of the sizes and proportions of animals or animal parts, has traditionally been considered a key element for breed characterization [1]. Such a determinant role not only relies on the implication of body conformation and dimensions in the definition of breed standards, but also on the functional classification of individuals depending on their better suitability for the development of certain tasks [2]. In these regards, zoometric analysis may help to detect differences among livestock populations or breeds which may be the source for niche specialization or exploration opportunities. As a direct consequence, zoometrics has been reported as one of the driving agents for the conservation (seeking a particular breed standard) and selection strategies (intraherd and interherd breeding criteria definition) that are eventually implemented [3] in local breeds.
Animal body measurements have been traditionally obtained manually with the use of a diverse range of instruments [4]. However, these tasks are not exempt of risks and inconveniences for both the animals and the workers, such as the increased stress induced in the animals or the errors in measurements due to the difficulty of maintaining the animals in a complete stationary position, among others. If sedation or the use of anesthetics is required for animal immobilization due to temperament issues, therefore becoming potentially hazardous to the operators, not only can the animal welfare be compromised, but this also makes the time and costs of the zoometric tasks increase [5].
For these reasons, several methods aiming at automatizing zoometric measurements collection in non-invasive, contactless, cheaper, and faster ways have been attempted in different species over the past decade [6]. Literature has reported the success of bidimensional and tridimensional image-based zoometric analysis methods in domestic and wild animal populations [7][8][9][10][11][12][13][14][15][16][17][18][19][20]. Among them, software-assisted digital imaging zoometrics offers a sound and solid alternative, which has not only solved the aforementioned issues, but also overcome the accuracy of traditional zoometric practices [5,6]. The evaluation of static images may reduce the biases derived from human data collection due to the animals' spontaneous movements. However, the use of high-resolution images may be needed, and inferring certain tridimensional measurements such as perimeters from bidimensional photographs is still critical. In this sense, mathematic modelling can help correct potential computational biases from both linear and tridimensional measures, while at the same time providing an economically affordable opportunity to perform zoometric comprehensive analyses, and the stress induced in animals at the time of being held and measured is substantially minimized [5].
Some minor species such as camels, despite being increasingly prevalent at livestock scenarios for their contemporary recognition as a sustainable species [21], have only anecdotally been attempted for zoometric evaluation and breed characterization [3]. Dromedaries or one-humped camels (Camelus dromedarius) are a typical element in the scene of developing economies and constitute the vast majority of the world's camel census. Due to the economical context in which they are normally evaluated and the lack of attention paid to them in the past [21], the scarce initiatives towards morphometrics phenotypic variability collection have mainly been implemented via on-field sampling [1,[22][23][24][25]. However, these practices are not free from challenges. Indeed, the combination of camels' large size and often strong temper may compromise the integrity of operators [3] and turn zoometric analysis into a dangerous time and a demanding human resources practice [16]. As a consequence, certain measurements, or medial regions like udders [26] and genitals [27], are routinely almost never registered due to the difficulty or danger that their access implies.
Contextually, for the safe accomplishment of on-field zoometric collection, camels need to be properly restrained for their secure handling [4], and the use of a wide diversity of measuring tools is compulsory. Hence, camels still lack contactless methods which may help to safely improve the efficiency and accuracy of zoometrics for camel breed characterization and the definition of adapted selection criteria for the maintenance of camel global genetic diversity [28]. This becomes even more crucial when the camel breed being considered is at risk of extinction and has a very defined applicability for tourism, work development, or therapeutic kinetics [29], as occurs for the Canarian camel.
The early development of image-analysis methodologies for measuring zoometrics in dromedaries is evidenced by the first body zoometric digital reconstruction [30]. As a result, the 3D modeling method is confirmed to be used as a remote method to extract morphological features of camels in a reliable manner. However, the experimental nature of this initiative and the time and costs needed for it to be implemented at a large scale make it compulsory to explore other alternatives. In this sense, a recent investigation has approached the accuracy of image analysis for the extraction of some linear zoometric measurements in dromedaries [31]. Nevertheless, the evaluation of tridimensional measurements (such as perimeters) is still misconsidered. Thus, the relatively scarce data that are obtained continue to be incomplete, given that relevant functionally important traits might not be recorded [32]. In this framework, the present research aims to develop a standardized and validated method for the comprehensive collection of both linear and tridimensional zoometric measurements in live dromedary camels by using 2D images. The combination of bidimensional digital imaging and logarithmically adjusted Bayesian regression methods provides a timely response to the need for the accurate registration of linear and tridimensional zoometric measurements. The method proposed not only offers a time and money affordable precise alternative that can be used as the main source for breed characterization and functional evaluation at a large scale, but also in a comparative manner which may be translatable to other camel populations.

Zoometric Parameter Definition
The bibliography on the topic was reviewed to obtain a comprehensive database of zoometric measurements in camels during the whole month of September 2019. Bibliographic analysis was performed using Google Scholar search engine (https://scholar.google.com/) (accessed on 1 September 2019), as suggested by other papers in which document library data extraction has been performed, due to the possibilities this search engine offers in regards to data extraction [21]. After this document search, six papers dealing with camel zoometrics, regardless of the measuring method used, published from 1994 to 2019, were found [33]. The list of measurements included in the aforementioned documents was completed with other variables relevant for camel functional development [3]. After a variable list was completed, a total of 30 zoometric measures were determined to be collected on-field and later extracted from digital images. Table 1 presents a description of the aforementioned zoometric variables.

Thorax and Dorsum
Chest width Distance between the medial point of the front of the forelimbs, measured at the base of their insertion into the torso.
Heart girth Circular perimeter of the chest measured directly behind the sternal callosity and before the hump. Height at withers (stature) Distance from the withers to the ground. Body length Distance between the shoulder to the point of the hip.

Hump
Hump-to-tail distance Distance between the most caudal point of the base of the hump to the base of the tail. Hump length Distance between the front and the back of the hump passing through the top of it.
Hump width Distance between the front and the back of the hump passing through one of its laterals.
Hump height Distance between the middle points of the base of the hump at each lateral, passing through the top of the hump. Hump girth Circular perimeter of the hump measured at its base

Rump and Tail
Rump length Distance from the coxal to the ischial tuberosity.

Rump width
Distance between the right and left coxal tuberosity. Tail length Distance from the base to the tip of the tail, excluding the tail skirt. Width at the base of the tail Distance between the most lateral points of the base of the tail.

Animal Sample
Zoometric record collection took place between September 2019 and August 2020 for 130 Canarian camel breed individuals (58 females and 72 males). Camels were located in three representative emplacements where Canary camels are bred: (Doñana National Park) Huelva (36.972330, −6.427498), Almería (36.902180, −2.429520), and Fuerteventura (28.186777, −14.158361) in Spain. Animals were clinically exanimated to ensure the proper condition of the animals, which enabled their participation in the study. Furthermore, to prevent bias issues derived from the sexual status of the animals at the moment of sampling, sexually immature individuals were discarded (below 3 years of age). Parallelly, only non-gravid she-camels were included in this study, since pregnancy may be a source of bias in zoometric measurements in the thoracoabdominal region [34]. Age or live weight did not normally distribute (p < 0.05). Live weight was calculated using the following formula (Equation (1)) by Boujenane [35].
Live Weight= 6.46 × 10 −7 (HW + ChG + HG) 3.17 (1) where HW is height at the withers, ChG is chest girth, and hg is hump girth, respectively. The sampling of each of the 30 zoometric measures taken from each animal was collected from its left side following the premises in Iglesias et al. [33] and Alhajeri et al. [3]. Females and males' age and live body weight descriptive statistics are represented in Figure 1.

Sampling
The end of the molting season, which is a six-to-eight-week period starting in late spring [36], was chosen as the sampling moment to prevent the bias which may potentially be ascribed to hair length and texture [8]. Figure 2 presents a flowchart summarizing the research methodology.

On-Field Zoometrics
Live animal-based measurements sampling took place with the animals holding a static upright position with their head naturally raised and in correct aplomb (parallel fore and hind legs perpendicular to the ground with lined toes). Animals were measured on a flat and hard ground surface. Measure collection was performed using a non-elastic measuring tape. All operators were trained. The first operator performed on-field measurement collection and digital image measurement extraction. The first operator was assisted by a second operator in zoometric measurements collection, while a third operator annotated the outputs of zoometric evaluation and held a one-meter measuring bar to be used as a reference for calibration on digital zoometrics extraction ( Figure 3). aplomb at a completely static position while keeping the photographic and camel midsagittal planes parallel; second operator was responsible to take photographs for digital imaging analysis; Second operator (Blue): responsible to take photographs for digital imaging analysis, and third operator (Orange): responsible to hold a one-meter measuring bar (yellow bar) near the animal body to be used as a reference for image calibration on digital zoometrics. During (B) on-field collection of zoometric data, first (Grey) and second (Blue) operators perform on-field measurement collection with a non-elastic tape (green tape) while third operator (Orange) annotates the outputs of zoometric evaluation. During (C) digital imaging analysis, first operator (Grey) performs digital image measurement extraction.

Digital Imaging
Three photographs (front, lateral, and back perpendicular to the camera) were taken per animal right before on-field zoometric evaluation. The second operator took these three photographs for digital imaging analysis (front, lateral, and back views). The third operator was in charge of holding a one-meter measuring bar at the same midline of the body to be used for reference calibration of distances on the computer measurement software for digital imaging zoometrics while taking the aforementioned photographs. The obtained images were digitally processed using Kinovea 0.95 (Free Software Foundation, Inc., Boston, MA). Zoometric linear measurements were obtained in pixels by drawing a straight line between two points in the picture and automatically converted into cm after calibration of the software using the measuring bar as a reference using the Calibrate option of the Line tool of the software [37]. Puig-Diví et al. [38] reported Kinovea software to be a valid and reliable tool which is able to measure accurately at distances up to 5 m from the object and at an angle range of 90 • -45 • .
Image collection was performed on an open hard ground and flat area. Light conditions were chosen so as for the animal not to be placed in a shaded area or in one in which light exposure may distort image capture. The animal color was considered to ensure background color did not lead to any measure distortion or measure misregistration. The camera was positioned at a standardized height of 1 m on a camera stand 4 m away from the camel center of balance. The aforementioned distance and height permitted framing of the animals being evaluated on the whole. We followed the premises in Iglesias et al. [33] to determine the proper aplomb of the animals and tracing standard lines on the ground before photograph taking to ascertain the animal was in the right position. Image capture was performed using a digital camera (Sony DSC-RX100 SENSOR CMOS Exmor 1.0 of 20.1 MP, F1.8-4.9, Zoom 20-100, Optical Zoom 3.6×, 3 LCD Screen Image stabilizer) in standard mode. Joint Photographic Experts Group (JPEG) compression format was used. One trained operator performed zoometric measurement digital extraction from photographs manually. Intraclass correlation coefficient (ICC), based on multiple paired Cohen's κ tests, was run to compare zoometric on-field analysis to digital imaging zoometric analysis. As suggested in Bunting et al. [39], the intraclass correlation coefficient (ICC) is a reference method to determine the reproducibility and reliability of numeric measurements organized into groups beyond a simple pairing, for example, different operators measuring the same variable in different animals or the same operator using different methods on different animals. In this study, we issued the equations, and equations were solved. Then, we used ICC to compare the results from model solving and real measurements to test for the reproducibility and reliability of models.

Scale Reliability and Repeatability: Cronbach's Alpha (Cα)
Instrument scale internal consistency was measured using Cα. Internal consistency, when applied to instrument comparison, is an estimate of 'reliability based on the average correlation among items within a test', with each of these items being each of the measuring methods, zoometric on-field evaluation, and zoometric digital imaging, in our case [42], and examines the degree to which such instruments measure the same characteristics or domains of knowledge [43]. Typically, internal consistency is measured by the calculation of a reliability coefficient [43], such as Cα.
In this context, Cα represents the reliability level of the instrument being compared to a reference instrument (zoometric digital imaging to zoometric on-field analysis as a reference) [44]. As a general criterion, George and Mallery [45] suggest the following recommendations for evaluating Cα coefficients: >0.9 is excellent, >0.8 is good, >0.7 is acceptable, >0.6 is questionable, >0.5 is poor, and < 0.5 is unacceptable. However, when comparing internal consistency between instruments, Pallant [46] reported that Cα value above 0.6 is considered a highly reliable and acceptable index [42]. Furthermore, retaining variables with values over 0.5 has been suggested due to their ability to explain data variability [47].
As suggested by González Ariza et al. [47], single measures of ICC determine how a single observation taken at random may correlate to another single observation, that is, in our case, how a zoometric measure from on-field evaluation may correlate with its paired counterpart from digital imaging. By contrast, average ICC and Cα determine how consistent the set of instruments being compared are on average. Consequently, in instrument comparison, average measures somehow prevent potential measuring errors affecting particular measurements, and as a result, report erroneously decreased instrument reliability and accuracy values for the instrument being tested (digital imaging zoometrics in our case).

Parametric Assumptions Testing and Approach Decision
The statistical approach was decided after parametric assumptions testing. The Shapiro-Francia W' test (for 50 < n < 2500 samples), Shapiro-Wilk test (for n < 50 samples), and Levene's test were used to determine whether normality and homoscedasticity parametric assumptions were met. The Shapiro-Francia W' test was run using the Shapiro-Francia normality routine of the test and distribution graphics package of the Stata Version 15.0 software (StataCorp, College Station, TX, USA) [48]. Homoscedasticity was run using the explore procedure of the descriptive statistics package in SPSS Statistics (Version 25.0, IBM Corp., Armonk, NY, USA) [49].

Perimeters and Girths Calculation
As suggested by Singaraju et al. [50] in their study performing ellipsoid biometric computations, Ramanujan's equation of ellipse model for ellipsoid perimeter (P) was used to fit zoometric perimeters and girths (neck girth (cranial, middle, and caudal thirds), heart girth, hump girth, thigh perimeter, hock perimeter, fore cannon bone perimeter, rear cannon bone perimeter, and foot perimeter, respectively) as follows; where a is the semimajor axis, b is the semiminor axis, and h is computed as follows h = (a − b) 2 /(a + b) 2 . This approximation is within about 5% of the true value, as long as a is not more than three times longer than b. Figure 4 schematically represents semimajor and semiminor axis digital imaging collection references for perimeters and girths computation. Ramanujan's equation of the ellipse model was applied using Microsoft Office Excel 2016 as suggested in the literature [51,52].

Bayesian Linear Regression Modelling and Natural Logarithmic Correction for Perimeters and Girths
Limited sample sizes derived from endangered populations may distort the distribution of variables that are presumably sampled from normally distributed populations, such as zoometrics. Such distortion may be ascribed to highly skewed data appearing as a result of valid outliers. A valid outlier would be an animal that is considerably smaller or larger than the rest, but which may still fit in the context of an endangered breed standard, even if it statistically distorts sample distribution properties.
If valid outliers are not present, skewness may not be forced towards one of the distribution ends, and applying the logarithmic bias correction may not be necessary. However, body score condition and live body weight have been reported to act as sources of bias when performing zoometric analyses in dromedary camels [35]. The bias effect of other factors, such as age or coat color, was also considered.
In such cases, regressing the variable measured in the field against the potentially originating factor of the bias during the collection of measurements from digital sources (live body weight, age, or coat color in our case) and the same variable measured with digital methods (given our aim is to perform a comparison between methods) may have to be performed to ensure the replicability between methods. This ensures that the accuracy of estimation of real measurements after digital measurements is maintained.
The effects derived from outlier distortion become stronger in perimeters, girths, or circumferences calculations given that these cannot be extracted from photographs straight away. The inability to account for such a bias while evaluating tridimensional zoometric parameters from bidimensional images makes it compulsory to apply correction methods. According to the Energy Information Administration of the United States [53], to correct for skewness derived bias in estimations, a bias correction should be issued. To solve this issue, the Energy Information Administration of the United States [53] suggests that sample size limitations could be buffered by regressing the log transformation of the variable to which regression models were initially aimed.
Once the variable measured in the field has been regressed against the potentially originating factors of the bias, and the same variable is measured with digital methods, bias correction may be obtained after summing the outcomes of this regression to the natural logarithm of the zoometric variable measured using digital methods.
In these contexts, the application of regression equations may be difficult due to the alteration of sample properties and ordinary least squares regression assumptions. However, these distribution distortions events may be saved using statistical alternative methods such as Bayesian linear regression, which are less sensitive to outliers and distribution alterations [54], and their estimations may be subject to wide confidence intervals.
Valid outliers' detection in our data sample was performed using the identify outliers procedure of the Analyze/Built-in analysis of the column analyses package of GraphPad Prism version 8.3.0. The ROUT method was applied to prevent the effects of outliers. The ROUT method combines robust regression and outlier removal and is based on the false discovery rate (FDR). A maximum desired FDR must be predefined (Q coefficient). ROUT method assumes all data except for outliers to be sampled from a Gaussian distribution. When data does not meet the aforementioned assumption, outliers may follow the same distribution as data.
The ROUT method strength to detect outliers was determined using the Q coefficient. Higher levels of Q are indicative of lower threshold strictness for outliers' detection, thus, the higher the outlier detection power, but also the higher the probability for a false outlier to be identified as a true one. Lower Q values set stricter thresholds for outlier definition, which consequently translates into lower power to detect real outliers, but also lower chances for false outlier consideration. A Q coefficient of 1% is the recommended threshold to be used as a default [55], given it implies a lower than 1% false discovery rate for outlier detection.
Consequently, after considering the presence of outliers, when evidence of a lack of acceptable fit (Cα < 0.600) between perimeter or girth on-field zoometric evaluation and digital imaging mathematical computation using Ramanujan's equation for ellipsoid perimeter was found, we used Bayesian linear regression modelling to correct measurement as a function of live body weight.
The presumably large dependence of zoometric parameters on age may have suggested the inclusion of age in the regression models, as according to Carlin [56], Bayesian inferences are sensitive to the dependence of variables on time (conditional on θ and x). In our case, zoometric analyses were performed when animals had reached the adult stage to prevent age-derived biases from occurring, as suggested by the lack of pieces of evidence or a significant effect of age (p > 0.05). Contextually, the age variation coefficient (CV) was 0.598. As a rule of thumb, a CV ≥ 1 indicates a relatively high variation, while a CV < 1 can be considered low, hence, the lack of a significant interindividual variability prevented age from being considered as a covariate.
Pieces of evidence for a non-significant bias effect of coat color effect were also reported (p > 0.05 and CV of 0.376). This lack of bias ascribed to coat color variability may derive from the fact that photographs were taken once the hair molting season had passed and an appropriate lighting scenario that remarked the topographic shadows of bone accidents (references for measurement) had been considered. Furthermore, animals were photographed and measured by the same people trained for this purpose, with the same technological equipment, and information regarding their coat color was registered to minimize the effects of the aforementioned potentially biasing factors.
Those zoometric perimeters and girths did not reach acceptable levels of scale reliability (Cα < 0.600), and hence were considered the dependent variables of Bayesian regression models to apply body live weight correction.
As suggested in Koehrsen [57], Bayesian linear regression uses probability distributions rather than point estimates. This means response, y, is not estimated as a single value, but is assumed to be drawn from a probability distribution. The output, y, is generated from a normal (Gaussian) distribution characterized by a mean (the transpose of the weight matrix multiplied by the predictor matrix) and variance (the square of the standard deviation σ, multiplied by the identity matrix, given it is a multi-dimensional model formulation). Hyperparameters' means and variance were used to obtain the best values of the hyperparameters of the prior distribution, as suggested by Kundu [58].
The objective of Bayesian linear regression is to determine the posterior distribution for the model parameters, rather than the best value for model parameters. Not only is the response generated from a probability distribution, but the model parameters presumably come from the same distribution. The posterior probability of the model parameters is conditional on the training inputs and outputs. In contrast to frequentist Ordinary Least Squares Regression (OLS), there is a posterior distribution for the model parameters proportional to the likelihood of the data multiplied by the prior probability of the parameters.
This implies two primary benefits of Bayesian linear regression: priors and posteriors. When there is knowledge, or a guess for what the model parameters should be, these priors can be included in the model (for example, the influence of live body weight on zoometrics). This approach contrasts the frequentist approach, which assumes everything there is to know about the parameters comes from the data. Indeed, in Bayesian regression, when there is no prior information known, non-informative priors for the parameters such as a normal distribution can be used.
Afterwards, posteriors, or the results of performing Bayesian linear regression, are a distribution of possible model parameters based on the data and the priors. Posteriors enable quantifying uncertainty about the model. Hence, the fewer data points, the greater the dispersion of posterior distribution will be. As the amount of data points increases, the likelihood washes out the prior, and in the case of infinite data, the outputs for the parameters converge to the values obtained from OLS.
To summarize, in Bayesian inference for linear regression, we use priors as initial estimates, and as we gather more evidence, testing our model against data (posteriors), the model supports or disproves our prior hypotheses. In practice, the evaluation of the posterior distribution for the model parameters is intractable for continuous variables when we implement Bayesian linear regression. Thus, sampling methods are used to draw random samples from the posterior distribution to approximate the posterior distribution to which it should be using Monte Carlo algorithms method and its variants.
The (Metropolis-Hastings) random walk algorithm, which uses Markov Chains to perform Monte Carlo estimate via the Gibbs Sampler algorithm, was used as aforementioned given a different prior to the default uniform prior specified in IBM SPSS Statistics Algorithms v. 25.0. by IBM Corp. [59] was selected. The random walk Metropolis algorithm is the preferable option for data imputation from the collection of Markov Chain Mote Carlo (MCMC), as suggested in MacKay [60], given that neither admissibility nor stability were selected.
The following general equation was used for each of the regression models defined in this study y i = X 1 β 1 + . . . X i β i + ε i , where i = 1, 2, . . . i is the ith number of factors; y i is the vector of records for the aforementioned dependent variables with dimension n (a total of 910 records, one record per each of the seven circumferences/perimeters which did not reach acceptable reliability levels and each of the 130 dromedary camels measured); X i is the appropriate incidence matrix for factors; β i is the standardized regression coefficients for the ith number of factors and covariates considered, respectively. The general regression equation used for each perimeter or girth was Y = Intercept + β on field perimeter/girth (cm) ·onfield perimeter/girth (cm) + β Live body weight (Kg) ·body live weight (Kg).
As Brewer [61] suggested, the intercept was necessary given that we used unstandardized coefficients. The magnitude of intercept confidence intervals was an empirical indicator of the need for its estimation. Residual effects (ε i ) normality was assumed as follows ε i X i N 0, σ 2 εi , where X εi is an identity matrix and σ 2 εi is the residual variance, re-spectively. Continuous predictor variable unstandardized coefficients were produced by the linear regression model using the independent variables measured in their original scales. As suggested by Hayes et al. [62], unstandardized coefficients (βi) can be defined as the average increase of βi units in Y associated with an increase of one unit in Xi maintaining the rest of the variables constant. Below, a detailed summary of the priors and posterior distributions used in this study is reported. The complete description of the algorithms used by SPSS to perform Bayesian inference on multiple linear regression models in this study can be found in IBM SPSS Statistics Algorithms v. 25.0. by IBM Corp. [59].
Quadratic approximation was discarded (even if it has been reported to be computationally faster in terms of discretization and computing the likelihood over all possible parameter combinations). Instead, the Markov Chain Monte Carlo (MCMC) approximation was used, as it does not assume the fact that the posterior distribution follows a normal distribution.
Bayesian linear regression analyses were performed using the linear regression package from the Bayesian statistics task of SPSS Statistics, Version 25.0, IBM Corp. [49]. The Bayesian linear regression test routine of the linear regression and related package of the Stata Version 16.0 software process was used to compute posterior distribution statistics for the factors considered. Afterwards, we evaluated the estimated effect of the factors considered in the resulting predictive models, its confidence intervals, and posterior distribution statistics to build linear regression equations, calculate digital perimeter/girth extrapolation, and ICC was run again. When Cα had not significantly improved over acceptable levels, the natural logarithm of each particular perimeter or girth was added to each aforementioned equation, and ICC was calculated again to ensure that reliability levels had been attained.

Jeffrey-Zellner-Siow (JZS) Mixture of g-Priors
The Jeffrey-Zellner-Siow mixture of g-priors [63] was used given it successfully satisfies several theoretical requirements such as the equality constraint on the test-relevant parameters, for instance of β, which leads to the null hypothesis H 0 = β = β 0 [64], as suggested by Heck [65]. Rouder et al. [66] and Liang et al. [63] also acknowledged the benefits of JSZ prior distribution. Contextually, conditional on the residual variance (σ 2 εi ), the JZS prior defines a multivariate Cauchy distribution for the slope parameters of the full model, as follows.
(β i σ 2 εi ) ∼ MVC 0 P , γ i 2 σ 2 εi C i −1 , which is defined by a P-dimensional zero vector (location vector) and a scale matrix. The constant γ i measures the magnitude of scaling and is a priori chosen by the user, the residual variance σ 2 εi , and the matrix C i = X i X i /N i , which is the covariance matrix of the centred design matrix X i .
JZS prior [66] is especially appropriate in Bayesian linear regression analyses given that it is symmetric and centered at zero, as explained by Bayarri et al. [67]. This means that positive and negative values of the parameters of the slope are a priori equally likely to occur. Moreover, JZS prior does not depend on the scale of the variables, factors, or covariates considered. Hence, the Bayes factor is scale-invariant, and outputs remain the same when the variables expressed in different units are evaluated, which is likely to occur in multifactorial studies.
Scaling the multivariate Cauchy distribution by the residual variance σ 2 εi ensures the achievement of such independence from the measurements of model elements (a priori, a larger residual variance implies larger slopes) and by the inverse of the covariance matrix C i (a priori, a covariate with a larger variance implies smaller slopes). In this context, the process of definition of scaled priors for unstandardized coefficients (β i ) equals the process of definition of priors for standardized coefficients (β * i ) [66]. Third, the scale parameter γ is fixed to a constant, hence prior beliefs are specified about the expected effect size remain constant as well. The IBM Corp. [59] algorithm guide, in its section for the algorithm of JZS prior for linear regression analyses, sets the default value of γ = 2 √ π = 3.5 to compute Bayes Factor. This reflects a belief of a priori medium effect size, which, for a single covariate x, implies a priori probability for the standardized regression slope β * i = β i · SD(x i )/σ i of 53.2% of being in the range (−0.50, +0.50). Rouder and Morey [68] explained other benefits from the choice of the JZS prior, for instance, its model selection aimed at consistency (this is that Bayes factor, goes to infinity as the number of observations N increases without bound-favoring the data-generating model) or consistency in information (the Bayes factor for a certain effect goes to infinity as the proportion of explained variance or R Squared (R 2 ) increases to 1). Additionally, Bayes factors for JZS prior are highly precise and relatively easy to compute [69], resulting in its wide applicability for the default t-test [70], ANOVA [66], and linear regression [65].

Bayesian Modelling of Factor and Covariate Effects (FCEBM)
Being y i , any of the effects of any of the independent variables (covariates) considered in this study (live body weight, age, and coat color), the posterior distribution of y i in the context of the data, D, is This means each model's average of posterior distributions is weighted by their posterior model probabilities. In the aforenoted equation, the posterior predictive distribution of y i given a particular model M i is, and the posterior probability of the model M i is given by where, is the integrated probability of the model M i , β i is the vector of parameters of the model M i , p(β i |M i ) is the prior density of β i under model M i , p(D|β i |M i ) is the probability, and p(M i ) is the prior likelihood that M i is the true model. The number of models (K) for a problem with P potential covariates can be enormous (K = 2 P in the absence of other constraints). Only a small number of these K models will be sufficiently supported by the data, and hence selected by SPSS for each of the P covariates. Gibbs sampling algorithm was used to estimate marginal posterior distributions of all unknowns.

Factors and Covariate Effect Bayesian Interpretation (CEBI)
The detection of issues before model estimation was evaluated using the checklist proposed by Depaoli and Van de Schoot [71]. Among the issues checked, we tested for those occurring after model estimation before interpreting results, in priors' influence comprehension, and after interpreting results for conclusion drawing. Interpreting the effect of each particular covariate (independent variables used in this study) was made as follows.
First, the posterior probability p[β * i = 0/D] expresses the probability that every single independent factor or covariate affects each particular dependent variable. Standard rules of thumb [72] for posterior probability interpretation are as follows: <50% evidence against the effect; 50-75% weak evidence; 75-95% positive evidence; 95-99% strong evidence; >99% very strong evidence, which is comparable to commonly used thresholds that define the level of significance of evidence using Bayes factor (BF) ( Supplementary Table S1).
Second, posterior distribution means determines the magnitude of the effect of every single factor and covariate. For metric covariates (continuous predictors) or the numeric variables used in this study, regression coefficients define the difference in the predicted value of the response variable for each one-unit change in the predictor variable, with all other predictors being constant. When dependent variables are metric, β regression coefficients are a measure of effect sizes by themselves.
Third, a 95% credibility interval suggests that a 95% likelihood for these regression coefficients (every single covariate and factor posterior distribution means) lies within the corresponding credibility intervals. A significant effect is reported when 0 is not contained within the credibility interval for each particular factor. In the present study, only live body weight significantly influenced zoometric measurements (p < 0.05).

Convergence Criterion
Iteration rounds continued until a tolerance convergence criterion of 10 −8 was reached [73]. After this, initial parameters were defined, and model fitting properties were analyzed. The maximum number of iteration rounds was 2000 for each analysis as stated in IBM SPSS Statistics Algorithms version 25.0 by IBM Corp. [59]. Such a convergence criterion was defined given its wide application in Bayesian ANOVA and linear regression analyses in limited sample sizes research contexts [74].

Model Validity and Explanatory Power of Present Data, and Predictive Power of Future Data
Validation and comparison of Bayesian models were described in Geweke [75]. Contextually, other authors [76] suggest model validation should base on models' mean square error (MSE). Although mean square residual or error (MSE) and minimum mean-square residual or error (MMSE) have been widely used to measure the closeness between a regression line and a set of points (model fit to explain data), mean square prediction error, or MSPE (=RSS/no. of observations), was used to measure error variation due to the MSE being influenced by the number of predictors [77] in reduced sample sizes cases [78,79].
The residual sum of squares (RSS) quantifies the amount of variability in a data set not explained by a regression model. That is, the RSS measures the amount of error remaining between the regression function and the data set, thus it essentially defines the ability of a certain regression model to explain or represent the data. Smaller values of RSS are indicative of better suitability of the regression function to model for the data it intends to model. Specific to Bayesian inference, Monte Carlo standard error (MCSE), which is defined as the standard deviation of the chains divided by their effective sample size, measures the accuracy of the chains. MCSE is the non-parametric or Bayesian counterpart of MSPE, and should be used as the validation criteria in Bayesian linear regression model comparison studies [80].
Bayes factor (BF) is an indirect measure of models' explanatory power to describe observed data. Larger values of BFs evidence higher probabilities for the combinations of the factors modelled to explain dependent variables. Supplementary Table S1 reports common thresholds to define the significance of evidence as suggested by Jeffreys [81] and Lee and Wagenmakers [82]. Bayesian R 2 is related to BF and can be considered as a data-based estimate of the fraction of variance explained for data. Parallelly, acceptance rate, efficiency, and Monte Carlo standard error (MCSE) were used to determine Bayesian methods' validity. Supplementary Table S2 presents the description and interpretation of model validity parameters. The predictive accuracy of the model's Bayesian statistics [83] can be estimated through posterior predictive checking [84].
Afterward, the Bayesian information criterion (BIC) or Schwarz information criterion (also SIC, SBC, SBIC) was calculated to determine the model predictive ability of new data as follows: where MCSE is the Monte Carlo standard error, N is the number of observations or records, and K is the number of independent parameters of the model. BIC was used to compare predictive power across models as it considers the statistical goodness of fit and the number of parameters to be compulsorily estimated to reach such fitness degree, as it penalizes every time the number of parameters considered increases [85,86]. As a result, BIC quantifies the balance between model fit and model complexity [87]. Lower BIC values mean that a particular model is a comparatively better predictive model than the rest. Contrastingly, Bayesian R 2 estimates the explanatory power of observed data.
In consequence, although the addition of "noise" variables to the fit, variables that explain small redundant fractions of variance, will slightly increase R 2 values, model predictive power will decrease (as denoted by its higher BICs). Indeed, as more variables are added to the model, its predictive accuracy decreases. This is, higher R 2 will also translate into higher-and therefore worse-BIC values.

Parametric Assumptions Testing and Approach Decision
All variables did not grossly meet the normality assumption (p < 0.01), respectively. Homoscedasticity was violated as well (p < 0.01), hence, a parametric approach was discarded, and Bayesian methods resulted in being the most preferable option. No likely outlier was detected, and therefore we preserved all observations for further analyses.

Initial/First Round of Method Comparison/Repeatability and Scale Reliability
Among the zoometric perimeters and girths measured, neck girth at cranial and middle thirds and hump girth reached acceptable levels of scale reliability (Average/Cα ≥ 0.600, Table 1). Therefore, their direct extrapolation after the application of Ramanujan's equation of the ellipse model for ellipsoid perimeter calculation was feasible and no further analysis was performed. The rest of the variables, which did not reach acceptable reliability levels, continued for further analyses (Average/Cα < 0.600, Table 2).  a First: Variables for which direct translation between methods was feasible; Second: Variables for which Bayesian regression modelling correction as a function of digital measurements and live body weight was performed (except for foot perimeter, for which live body weight did not report a significant effect (p < 0.05)); Third: Variables for which Bayesian regression modelling correction as a function of live body weight and natural logarithm correction addition was performed.

Bayesian Linear Regression Modelling and Second Round of Method Comparison/Repeatability and Scale Reliability
After the initial round of method comparison/repeatability and scale reliability, neck girth at caudal third, heart girth, thigh perimeter, hock perimeter, fore cannon bone perimeter, rear cannon bone perimeter, and foot perimeter, respectively, did not reach acceptable levels of scale reliability (Average/Cα < 0.600, Table 2). Hence, we subjected them to Bayesian Regression Modelling Correction as a function of live body weight (except for foot perimeter, for which a nonsignificant effect of live body weight was reported p < 0.05), and the second round of ICC computation was performed. We did not consider additional factors such as coat color or age, given no evidence for their significant effect was found (p > 0.05). Table 3 summarizes Bayesian unstandardized linear (β) regression coefficients posterior distribution statistics for each of the variables considered. Bayesian determina-tion coefficients (R 2 ) or percentages of variance captured by each of the models and their respective BF are provided in Table 4. Models were considerably more probable than those which only comprised the intercept. Posterior predictive p values (pppvalues) for models were 0 < pppvalues < 1 and around 0.500, so therefore model fit was ensured (Table 5).

Natural Logarithmic Correction for Perimeters and Girths and Third Round of Method Comparison/Repeatability and Scale Reliability
Afterwards, all the variables except for thigh perimeters had not reached acceptable levels of scale reliability (Average/Cα < 0.600, Table 2). Therefore, natural logarithm was added to their equation, and the third round of ICC computation was used.

Final Outputs for Method Comparison/Repeatability and Scale Reliability
After Bayesian linear regression modeling and natural logarithmic correction, all methods proved to be highly repeatable and reliable, as suggested by the values of ICC and 95% CI for average measurements (Average/Cα) for all the zoometric measurements that were made. Finally, the resulting equations are reported in Table 4. Linear zoometrics always reported excellent repeatability values over 0.836 (for sole length average measurements ICC) with high 95%CI lower bounds over 0.768. Although perimeters and girths reported lower average measurements ICC, values were always higher than 0.654 (for foot perimeter average measurements ICC), which suggested good repeatability, while also presenting a reasonable 95%IC lower bound over 0.511. A summary of the results obtained after the three rounds of ICC calculation (single and average measurements/Cα) and 95%CI is reported in Table 1. Dispersion statistics for the measurements collected on-field and those extracted from digital imaging methods after correction are reported in Supplementary  Table S3.

Discussion
Apart from the time reduction, personnel requirements, and handlers' safety improvement implicit in digital imaging zoometrics, a lower expected random error is obtained when these are compared to on-field measurements. This improved accuracy, reliability, and repeatability is not only supported by our results (Tables 2 and 5, Supplementary  Table S3), but also by literature references highlighting a greater dispersion of the zoometric data collected with traditional instruments on-field against those extracted from digital image-based methods [5,8,88].
This may be ascribed to the difficulty for the animals to remain completely immobile during on-field measurement collection, which could add undesired noise to the data. As a drawback of image-based zoometric measurements, rigorous efforts may need to be made by operators to hold animals in the correct aplomb at a completely static position while keeping the photographic and camel midsagittal planes parallel [6].
In line with these preliminary arguments, collateral factors such as the image resolution, the accuracy of the distance measurement from the camera to the animal (focal plane), hair length and extension, coat color, and the fat condition of the animal, may need to be controlled before and during the zoometric measurements extraction from digital images to ensure that data depicts the real morphometry of the individuals as far as possible [8,19,[89][90][91].
The correspondence between on-field and image-based zoometrics methods was evaluated via the average measurements intraclass correlation coefficient (ICC)/Cronbach's alpha (Cα) [92]. Method correspondence was considered acceptable when a limit value of 0.6 ICC/Cα was reached [93]. This limit was attained for all the measurements except for 'Neck girth: caudal third', 'Thoracic girth', 'Thigh perimeter', 'Hock perimeter', 'Fore cannon bone perimeter', 'Rear cannon bone perimeter', and 'Foot perimeter' at the first round (direct extrapolation). These findings were compared with those of Pezzuolo et al. [94] and Pérez-Ruiz et al. [89], who reported a particularly large relative error in digital threedimensional zoometric measurements when compared to on-field three-dimensional measurements. Hence, the need for a correction methodology for digital three-dimensional measurements was suggested.
The basis for this lack of initial correspondence (in the first round) between threedimensional measurements may be based on the fact that for their on-field collection, the participation of at least two operators is generally required. This, in turn, results in an unbalanced placement of the corresponding measuring instrument at different transverse planes in the homologous contralateral regions. By contrast, linear measurements are generally made by the same operator who controls the exact position of the instrument used for the measurement and is only performed on one side of the animal, which reduces the probability of finding morphological variance throughout the measured area.
As a result, the second round of correction was necessary to enhance method correspondence for the measurements for which acceptable levels had not been attained in the first round. A Bayesian linear regression was applied between manual and non-contact measurements, considering body weight as the most critical influencing parameter.
Evidence of body weight is a source for such quantitative differences between onfield and digital imaging-based three-dimensional measurements found, as suggested by the lack of evidence of significance and CV values. According to literature, the circular perimeter of the area to be measured can intrinsically vary along its length depending on the regional adiposity [95] and the muscular development of the region [96], which are both directly related to body weight. Indeed, Boujenane [35] referenced that such a strong relationship between zoometry and weight could have been expected, since the estimation of the latter relies on chest and hump girth apart from height to withers via its estimation formula.
Once the digital measurement was multiplied by the coefficients obtained and the intercept value-added, the ICC between on-field and the corrected digital measurement improved only for the variable of 'Thigh perimeter'. Hence, the third round of correction was necessary to improve correspondence for the remaining measurements for which acceptable levels had not been reached. The correspondence between on-field and digital measurements was achieved after applying an additional correction consisting of the sum of the natural logarithm of the digital measurement (third round).
The fact that a linear regression model was enough to obtain a corrected digital measurement with acceptable correspondence for on-field thigh perimeter could be explained under the assumption that this area, given its location, hardly varies in its muscular tonus if no strong movements are carried out. This lack of variability does not distort skewness values; hence, the logarithmic correction does not improve correspondence values.
By contrast, for 'Neck girth: caudal third' and 'Thoracic girth', continuous respiratory movements may cause slight perimeter variations in the photographs taken during on-field evaluations, altering skewness properties and perhaps making it necessary to consider a natural logarithmic correction factor. This could also be applied to the rest of the variables, that is, 'Hock perimeter, 'Fore cannon bone perimeter', 'Rear cannon bone perimeter', and 'Foot perimeter'. The implicit difficulty to restrain camels in a completely static position may induce temporary size changes in these local areas. Specifically, the pressure that the fore and rear limbs have to cope with depends on the particular position of the animal aplomb at the moment of measure collection (photo taking).
This may derive from the fact that the musculoskeletal systems involved in maintaining the posture are incapable of maintaining limbs in a perfect and unique stationary position, with the consequent increase in local muscle temperature due to rising blood flow and the transient oscillation of the distal limbs to the resonant frequency of heartbeats producing muscle tremor [97,98]. This translates into a source of increased skewness, data dispersion, and noise, which is corrected, as suggested by literature, via the addition of natural logarithms of the zoometric measurement to the equation.

Conclusions
The efficiency, repeatability, and accuracy of the on-field zoometric methods that are routinely considered for zoometric characterization in dromedaries can be enhanced through the use of digital imaging techniques. On-field and digital imaging zoometric extrapolation results in comparatively improved performance of the latter, due to the reduced dispersion of the data extracted from photographs. However, mathematical correction methods may be needed for three-dimensional measurements to reach acceptable levels of on-field/digital imaging methods correspondence. Body weight is the main source of bias in perimeter digital zoometric extraction if other factors such as age and coat color are controlled. All zoometric perimeters collected on-field should be regressed against body weight, except for foot perimeter, to ensure acceptable levels of correspondence are attained. Additionally, we summed the natural logarithms of the particular measurements which could be affected by spontaneous movements (breathing or proprioception) to correct for derived skewness distortion. The present method offers a time and economically affordable alternative, which saves the drawbacks derived from the lack of consideration of three-dimensional measurements that routinely occur in digital imaging zoometrics. The high degree of correspondence between methods makes this tool valid for its standardized implementation in camel zoometric characterization, while it may be translatable to other species.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/math10193453/s1, Table S1. Commonly used thresholds to define significance of evidence through Bayes factor (BF); Table S2. Model validity and accuracy parameters definition and interpretation; Table S3. Descriptive statistics and dispersion measurements for digital imaging and on-field zoometric variables in Canarian camels.  Data Availability Statement: Data will be made available from corresponding authors upon reasonable request.