A Bayesian analysis of variables causally associated with hair cortisol concentration in dogs with obesity

Objective: To identify clinical variables causally associated with hair cortisol concentration (HCC) in dogs with obesity using a Bayesian analysis.

Study design: A retrospective analysis of clinical data and samples gathered from a cohort of dogs with obesity undergoing therapeutic weight reduction.

Methods: Hair was clipped from two sites (jugular groove, antebrachium), from dogs attending a specialist obesity care clinic, and combined before storage at −20 °C. Hair cortisol concentration was measured by liquid chromatography mass spectrometry. Causal associations between HCC and different clinical variables were assessed, informed by a directed acyclic graph. Variables assessed included age, sex, breed group, coat colour, body fat mass, weight reduction and the presence of comorbidities. Statistical analyses involved Bayesian multi-level modelling, with the magnitude of causal effects estimated using simulation from the posterior probability distributions.

Results: In total, 73 hair samples were collected from 52 dogs, with 31 providing single (before weight reduction) and 21 providing paired samples (before and after weight reduction). Dogs were of different ages, sexes and breeds, with most (44/52) having one or more comorbidities including orthopaedic, skin, cardiorespiratory, dental and neoplastic diseases. Mean HCC was 10.4 (standard deviation 19.52) pg/mg (logHCC 1.3, standard deviation 1.36). Bayesian multi-level models provided strong evidence that greater body fat percentage (98% probability) and presence of one or more comorbidities (>99% probability) were causality associated with increased HCC. Causal associations with other variables including, age, breed, sex, coat colour and season were less convincing.

Conclusion: Greater adiposity and having at least one comorbidity are causally associated with hypothalamic–pituitary–adrenal axis upregulation in dogs with obesity. Mechanisms warrant further investigation.

Introduction

Cortisol, the product of hypothalamic–pituitary–adrenal (HPA) axis activation, is a key hormone released in response to stress (1). Whilst acute cortisol concentrations can be measured in blood, saliva, or urine, such concentrations can be are influenced by circadian rhythms and acute stress at sampling (2, 3). Since cortisol is gradually incorporated into the keratinizing shaft during hair growth, hair cortisol concentrations (HCC) reflect its cumulative secretion over a period of weeks to months (1). Therefore, HCC has emerged as a non-invasive biomarker of HPA upregulation in both human and veterinary research (1, 2, 4). The method has been validated in dogs, showing good correlation between HCC and salivary cortisol concentration, and it has been applied to evaluate sustained cortisol output in various conditions (2, 4). For example, HCC is increased in dogs (5) and people (6) with spontaneous hyperadrenocorticism. Measurement of HCC is also used as a marker of possible stress in dogs experiencing different lifestyles or environments, such as working dogs and in dogs subjected to prolonged stressors in experimental settings (7). Therefore, HCC has the potential to capture chronic HPA axis upregulation stemming from both external and internal stressors. However, a standardised sampling protocol is required to ensure representative HCC measurements, controlling for possible confounding factors including coat colour, sampling season and region of the body from which hair is taken (8).

Chronic illness can also induce HPA upregulation in various ways including response to infection, chronic pain and due to altered immune function (9, 10). Obesity is a prevalent chronic disease in pet dogs, being associated with multiple comorbidities (11), a poorer quality of life (12), a shortened average lifespan (13) and also both functional and metabolic disturbances (14–17). Therefore, it is plausible that there may be HPA upregulation in dogs as a consequence of obesity. Indeed, in humans, there are positive associations between HCC and obesity metrics, including body weight, body-mass index and central fat distribution (18). In one weight maintenance trial, baseline HCC was associated with body mass in some cohorts, whilst increased HCC over 12 months predicted greater subsequent weight variability (19). One previous study has measured HCC in dogs with obesity, with little difference seen before and after a short weight reduction intervention involving lead walking (20); however, HCC tended to be greater in dogs undertaking the most exercise, possibly suggesting HPA upregulation due to physiological stress.

The term ‘causal inference’ is used to describe a process of analysing data to draw conclusions about causal relationships, and has most applicability when controlled experiments, including randomised control trials (RCT), are either impractical or unethical (21). Whereas traditional epidemiological methods focus on associations amongst variables, the aim of causal inference is to disentangle true causal effects from spurious correlation, including from confounding or bias. Although not yet fully exploited in canine research, causal inference has been used to assess risk factors for small intestinal dehiscence after surgery (22), unsuccessful dog ownership (23), early-onset incontinence (24) and patterns of physical activity (25).

Bayesian statistical analysis is an inferential framework based on Bayes’ theorem, in which prior knowledge or beliefs are updated after considering newly-observed data, enabling the probability of parameters or hypotheses of interest to be estimated (26). A crucial distinction from the frequentist approach is to consider the entire (posterior) probability distribution of an unknown quantity of interest, thereby better accommodating uncertainty in scientific research. Bayesian methods are particularly useful when sample sizes are small, when complex modelling is required (for example, hierarchical models) and when the degree of uncertainty needs to be quantified (27). Challenges with the computational complexity of Bayesian methods have recently been overcome by the development of Markov Chain Monte Carlo (MCMC) methods and availability of freely available online statistical software such as R (28); therefore, the use of Bayesian approaches is now feasible in biomedical research (29), including applications to veterinary species (30–34). Although most veterinary researchers are not familiar with Bayes theorem, it should arguably be intuitive to practising veterinarians because the method of Bayesian updating (sequentially updating initial beliefs as new evidence becomes available) is analogous to the process by which veterinarians investigate and make diagnoses in their patients (35). With this background in mind, the aim of the current study was to apply causal inference within a Bayesian workflow to investigate associations between body fat mass and other variables on HCC in dogs with obesity.

Methods

Animals

All participating animals were referred to a specialist obesity care clinic for dogs and cats (Royal Canin Weight Management Clinic, University of Liverpool, Neston, UK) for investigation and management of obesity or obesity-related disorders between May 2018 and February 2023, with all successful weight reduction interventions completed by September 2023. To be eligible, animals had to have had at least one adequate sample of hair available for cortisol analysis, and had to have reached an end point for their therapeutic weight reduction protocol, as described for similar studies in dogs (36); in this respect, some dogs completed their protocol and reached target weight, whilst others stopped prematurely, with the reasons recorded. Dogs diagnosed with either hypothyroidism or hyperadrenocorticism were not eligible given potential effects of these diseases on hair growth and HCC (2). However, having other comorbidities was not a reason for exclusion, and dogs were also not excluded because of the treatment they had received.

Collection of hair samples and measurement of hair cortisol

Hair samples from each dog were obtained by clipping two sites (jugular groove, antebrachium), whilst preparing for jugular venepuncture and venous catheterisation, respectively. Samples were collected before undertaking any clinical procedures including sedation for dual-energy X-ray absorptiometry (DXA). The quantity of hair collected depended on coat density, although this was always sufficient for HCC measurement. Hair from each site was combined, placed in an individual sealable plastic bag (Ziploc; S. C. Johnson, Wisconsin, USA) and then stored at −20 °C under dark conditions. All samples were subsequently shipped to a commercial laboratory (Dresden Lab Service, Dresden, Germany), with HCC measured in powdered samples by liquid chromatography mass spectrometry. This method has previously been used to measure HCC in dogs (37); the reported quantification limit was 0.09 pg/mg, whilst reported inter- and intra-day coefficients of variation were <10% (38).

Measurement of body weight and body fat percentage

Body weight was measured by electronic weigh scales, which were regularly calibrated using test weights (2–50 kg; guaranteed to be accurate to within ≤0.5%; Blake and Boughton Ltd., Thetford, UK). In most dogs (26/31 dogs providing single samples; 17/21 [before weight reduction] and 12/21 [after weight reduction] dogs providing paired samples), body fat percentage was analysed using fan-beam DXA (Lunar Prodigy Advance; GE Lunar; Madison, USA), calibrated on a weekly basis using a phantom supplied by the company, in conjunction with a bespoke computer software package (Encore 2004, 8.70.005; GE Lunar) (39). Dogs were either sedated (if DXA alone was performed) or anaesthetised if required for additional procedures, and scanned in dorsal recumbency, as described in a previous study (39).

Therapeutic weight reduction protocol

Full details of the weight reduction protocol used have been published in previous research (36, 40), although all dogs followed a partial weight reduction plan, meaning that the target weight set was deliberately greater than their ideal weight range, as described in a previous study (41). Briefly, at the first visit, patients were weighed, their body condition score (BCS) recorded and, in most dogs, body composition was also measured by DXA (see above). Health status was determined by routine haematology, serum biochemistry, free thyroxine measurement and urinalysis. If necessary, additional diagnostic investigations (e.g., diagnostic imaging, additional laboratory investigations) were performed to determine the status of any comorbidities. A tailored therapeutic weight reduction protocol was then formulated for each animal, again as described in a previous study (36, 40, 41). Briefly, animals were fed high protein, high fibre dry or moist therapeutic diets (Table 1; Royal Canin, Aimargues, France), with some dogs consuming the dry diet exclusively, and the remainder being fed a combination of wet and dry food, the choice of which depended on owner and animal preferences. The methods used to calculate initial food allocation have again been described in a previous study (41). In addition to advice about feeding the therapeutic diet, owners also received tailored advice on lifestyle alterations to assist the weight reduction process. This included a physical activity plan, tailored to owner circumstances, individual animal factors and the presence of comorbidities. Advice could include recommendations about play activity, walking, running, agility training and hydrotherapy.

Table 1

Table 1. Average composition of the therapeutic diets used for weight reduction in 298 dogs with obesity.

After the initial visit (V0, before therapeutic weight reduction), animals were reassessed every 7 to 21 days to have their body weight measurements taken, and changes were made to the dietary and exercise plan if necessary. In dogs that reached their target weight, a final evaluation was conducted (V1) after therapeutic weight reduction (median follow-up 313 days; range 116 to 1,609 days). Health status was determined based on physical examination, haematology, serum biochemical analysis and urinalysis. Body weight and body condition were recorded, and body composition was reassessed by DXA.

Determining sample size

We included as many dogs as possible from those seen during the collection period, with final numbers being equivalent to, or exceeding those used in many previous veterinary studies where HCC was measured (42–45). Arguably, a formal sample size calculation is less critical when using Bayesian methods (46); such analyses automatically account for any uncertainty arising from the sample size because the entire posterior density distribution is reported, with wider intervals reflecting greater uncertainty.

Causal model development

A variety of factors may induce HPA upregulation and, therefore, increased HCC, in pet dogs with obesity. These relationships were illustrated by constructing a directed acyclic graph (DAG), a graphical tool that provides a bridge between substantive theory (in the form of a scientific model) and statistical analysis (47). Such graphs can identify potential sources of bias enabling statisticians to decide which additional variables should be selected for adjustment in a particular statistical analysis (48, 49). The DAG developed in this study was created using online software [DAGitty software, version 3.1 (50)], and its design was informed by relevant canine and comparative literature. With this software, causal pathways (that directly or indirectly connect the causal variable to the outcome variable), confounding variables and confounding (a.k.a. ‘backdoor’) pathways within a DAG (51) can readily be identified. To obtain an accurate estimate of a causal association, all backdoor pathways must be closed, whilst not closing causal pathways; this is done by ‘conditioning on’ (i.e., including in the model) a set of a variables that are not descendants of the causal variable of interest; doing this should block all backdoor paths, thereby fulfilling the so-called ‘back door criterion’ (51).

The final DAG is shown in Figure 1, whilst graphical representations of all adjustment sets are provided in the Supplementary File 1. As can be seen, this approach enabled appropriate sets of adjustment variables to be identified for modelling causal associations for most predictor variables, e.g., age, sex, breed, coat colour, body fat, season of sampling and comorbidity. However, appropriate adjustment sets could not be identified for both the ‘successful weight loss’ (i.e., comparing dogs completing a period of therapeutic reduction with those not completing) and ‘visit’ (before vs. after therapeutic weight reduction) variables. This was because, in the causal model, there were backdoor pathways that could not be closed on account of unmeasured confounding (‘backdoor criterion’ not fulfilled). These unmeasured variables included owner or environmental factors affecting feeding and lifestyle behaviours that might plausibly affect body fat mass, weight loss success and attendance at follow-up visits (Figure 1; Supplementary File 1).

Figure 1

Figure 1. Directed-acyclic graph (DAG), created using the online resource: https://www.dagitty.net, and based on the scientific model which was used to inform the statistical analyses. The outcome variable (marked “I”) is hair cortisol, whilst other variables in the model are either observed (blue circles) or unobserved variables (grey), whilst arrows indicate a causal association between one variable and another. This DAG was used to determine adjustment sets for each of the final models, as shown in the Supplementary File 1. Possible unmeasured variables influencing the ‘weight loss success’, ‘visit’ and ‘body fat variables could include owner factors (e.g., owner attitudes and behaviours in implementing a therapeutic weight reduction protocol), environmental factors (e.g., the living environment of the dog) and possible impacts from the COVID-19 pandemic.

Data handling and statistical analysis

Statistical modelling strategy

A Bayesian workflow was chosen for statistical analysis which, as far as possible, has been conducted and reported in compliance with the Bayesian analysis reporting guidelines (52). There were several reasons for selecting such an approach. First, Bayesian methods are particularly well suited to the type of model used (multi-level linear models), and the flexibility they allow in specifying models that are appropriate for the data (52). Second, they perform well when causal associations are being estimated, involving statistical analyses informed by a pre-specified scientific model based on a DAG (49). Third, Bayesian analyses are computationally robust and are better able to handle uncertainty when making predictions in the face of small sample sizes (49).

A fourth advantage of a Bayesian approach is the requirement that all assumptions, both scientific and statistical, be clearly and openly stated in advance, ensuring that a scientist ‘shows their working’. Scientific assumptions include the components of the scientific model, as detailed in a DAG, outlining a scientist’s understanding of the data-generating process (see causal model development section). Statistical assumptions are those made when selecting the type of statistical model, deciding on an appropriate likelihood function and in choosing appropriate priors (see the Variables, likelihood function and model parameters and Selection of prior distributions sections below). Not only must such assumptions be stated and justified in advance, but they should also be considered when interpreting results, thereby reinforcing the notion that “posterior inferences are only as good as the model and experiment that produced the data” (26). A fifth advantage is the use of hypothesis testing in Bayesian inference; rather than being limited to null hypothesis significance testing (e.g., determining evidence against a null hypothesis) as with frequentist statistical methods, probability estimates directly in support of a particular hypothesis can be made. Such hypotheses are more flexible and more intuitive to readers. Finally, as emphasised in the results section, for a Bayesian analysis, it is the entire probability distribution that matters, not simply a point estimate (e.g., mean or median) or whether an arbitrary threshold for statistical significance (e.g., p < 0.05) was reached. Therefore, interpretations can be more nuanced, better reflecting the uncertainties of the scientific process and any study findings.

Dataset, variables assessed and missing data

The dataset on which all statistical analyses were conducted is provided in the online Supplementary File 2. Continuous baseline data are summarised as either mean [standard deviation (SD)] or median and range, whilst categorical data are reported as a number (percentage). Baseline variables recorded were age (in years), breed group [mixed breed (reference category), cavalier king Charles spaniel (CKCS), pug, retriever, other], coat colour [light (reference category), mixed, dark], sex [female (reference interval), male], comorbidities [no (reference category), yes], body fat percentage and season of sampling [spring (reference category), summer autumn, winter]. As discussed above, body fat percentage data were unavailable from 18 visits, including 9 initial visits (before weight reduction; 5 from dogs providing single samples; 4 from dogs providing paired samples) and 9 visits (after weight reduction; all from dogs providing paired samples). These data were considered ‘missing completely at random’ because the reason they were missing had nothing to do with observed and unobserved data (53). Otherwise, there were no missing data for any other variable assessed (Supplementary File 2).

Statistical software

Statistical analysis was performed using an online open-access statistical language and environment (R, version 4.4.3) (28). The additional packages used for data wrangling and visualisation were: ‘dlookr’ [version 0.6.3; (54)], ‘dplyr’ [version 1.1.4 (55)], ‘ggplot2’ [version 3.5.2; (56)], ‘psych’ [version 2.5.3; (57)], ‘reshape’ [version 0.8.9; (58)], ‘readxl’ [version 1.4.5; (59)], ‘tidyverse’ [version 2.0.0; (60)], ‘vioplot’ [version 0.5.1; (61)] and ‘vtable’ [version 1.4.8; (62)]. The additional packages used specifically for Bayesian modelling were: ‘brms’ [version 2.22.0; (62, 63)], ‘bayesplot’ [version 1.12.0; (63, 64)], ‘bayestestR’ [version 0.16.0; (65)], ‘loo’ [version 2.8.0; (66)], ‘priorsense’ [version 2.8.0; (67)], ‘rethinking’ [version 2.42; (68)] and ‘rstantools’ [version 2.4.0; (69)].

Variables, likelihood function and model parameters

The primary outcome variable was HCC, which was logarithmically transformed (logHCC) and standardised prior to analysis. Standardisation involves first centering the data (subtracting the mean) and then dividing by the SD, a process which was also applied to all continuous predictor variables. Such standardisation has several advantages; it ensures that the intercept corresponds to the mean value (set to 0), and also makes estimates of beta coefficients (a coefficient that quantifies how much the outcome variable changes for single-unit change in the predictor variable) easier to understand; by standardising, each beta coefficient then reflects the change in outcome variable for a 1-SD change in the predictor variable. As well as improving efficiency of the MCMC algorithm used in computation, it makes prior probabilities easier to set since positive and negative values represent positive and negative effects, respectively.

Separate models were constructed with logHCC as the outcome variable and the following causal predictor variables: age, sex, breed group, coat colour, season of sampling, presence of a comorbidity and body fat. The possibility of reverse causality, in the association between body fat and log HCC, was tested in a separate model; for this model, body fat as the outcome variable and log HCC as the causal predictor variable. For each model, adjustment variables were included as determined from the DAG (Figure 1; Table 2; Supplementary File 1).

Table 2

Table 2. Prior and likelihood specifications of all causal models used in analysis, along with the sets of adjustment variables included, as determined from the directional acyclic graph.

The statistical analyses used were multi-level Bayesian models, which included dog as a grouping variable (to account for multiple samples from some dogs). Such Bayesian models have the advantage of allowing greater flexibility in the desired model structure. Exact details of the variables and parameters used for each multi-level model are shown in Table 2, whilst full details of the statistical workflow (including code used, statistical outputs and graphs) are available online: https://github.com/AliG71/hair_cortisol.

Even after logarithmic transformation, data for the outcome variable (logHCC) remained modestly right-skewed. To determine the most appropriate likelihood function, preliminary models with different likelihood distributions (e.g., normal, Student’s t and skew-normal) were compared by leave-one-out (LOO) cross-validation using the ‘loo’ package (66). For each model, the leave-one-out information criterion (LOOIC) was calculated, which estimates how well a model can predict future data from the same distribution as the observed data, with smaller values indicating better predictive performance (66). A skew-normal distribution was ultimately chosen given its superior performance, as well as better ability to replicate the distribution of the logHCC data (Figure 2). This distribution is a generalisation of the normal distribution (with mu [mean] and sigma [standard deviation] parameters) but includes an additional ‘shape’ parameter (alpha) to allow for asymmetry, whereby positive and negative values indicate positive and negative skewness, respectively (70). Like the Student’s t distribution, it can be more robust to outliers, not least where there is asymmetry around the mean. All final models fitted the data distribution well, as shown by the posterior prediction checks outlined below (Supplementary File 3).

Figure 2

Figure 2. (a) Comparison of observed data [log hair cortisol concentration (log HCC)] versus simulations of normal [blue dotted line; mean 0, standard deviation (SD) 1], and skew-normal (red line, mean 0, SD 1, alpha 4) distributions. Note that the observed data have been standardised by subtracting the mean and dividing by the SD. The simulation of the skew-normal distribution better reflects the shape of the observed data. (b,c) Comparison of the empirical distribution of the observed outcome variable (logHCCC; thick dark blue curve) with the distributions of many replicated data sets (thin light blue curves) drawn from the posterior predictive distributions of body fat models that use either a normal (b) or skew-normal (c) likelihood function. Although the model with the normal likelihood function does a reasonable job of predicting the shape of the observed data, the fit is better when a skew-normal likelihood function is used. Specifically, the predictions of the left tail are better aligned (closer to parallel), with means of the predictions clustering around the mean of observed data, and the right-skew is better captured.

Selection of prior distributions

Prior distributions were selected for all parameters of each model, with the overall aim being to ensure they were weakly regularising, adjusted to ensure that that pre-data predictions would span the range of scientifically plausible outcomes. This was confirmed by graphical visualisations and prior predictive simulations from models that sampled from the prior probability distributions only (see below). Justification for the choice of each prior is provided in the Supplementary File 4, whilst details of the final prior distribution choices for each model are shown in Table 2. In most cases, neutral priors were chosen (including for body fat percentage), except for comorbidity and coat colour because previous scientific evidence suggested ‘informed priors’ to be more appropriate (Supplementary File 4). Further, a difference in the variance of HCC between dogs with and without comorbidities was evident (Figure 3) and, therefore, a prior for sigma (standard deviation for the likelihood function) was not included in models that included comorbidity, either as the predictor or an adjustment variable. Instead, sigma was estimated as a parameter within the model, using comorbidity as a single predictor variable (Table 2). Simulated prior probability distributions for the body fat and comorbidity beta coefficients are shown in Figure 4.

Figure 3

Figure 3. Combined box-and-whisker and dot plot of the observed log hair cortisol concentration (logHCC) data, stratified by the comorbidity variable. The logHCC data have been standardised by subtracting the mean and dividing by the standard deviation. The thick horizontal line represents the median, whilst the upper and lower limits of the box represent the inter-quartile range (IQR). The upper and lower whiskers extend as far as the largest or smallest, respectively, values that are no further than 1.5 × IQR from the IQR, whilst outlying points are shown as small black dots. Given the difference in variance between groups, it was necessary to build some models that accounted for unequal variance amongst dogs that differed by comorbidity status. Affected models included those testing causal associations between comorbidity and HCC and also body fat and HCC (since comorbidity was included in the adjustment set).

Figure 4

Figure 4. (a) Simulated prior probability density distributions (blue lines) for the beta coefficients for causal associations between body fat (a) or comorbidity (b) and log hair cortisol concentration (logHCC). The thick, red dotted line represents the mean of the density distribution, whilst the solid, thin grey line is intersects the x-axis at zero, indicating a value where the beta coefficient would be neither positive of negative. The body fat prior was assumed to be a normal distribution with a mean of 0 and a standard deviation (SD) of 0.5; the comorbidity prior was also assumed to be a normal distribution, but with a mean of 0.25 and a SD of 1. The impact that such prior probability distributions can have in regularising models is illustrated in Figure 8.

Model computations and diagnostic checks

Bayesian analyses were computed using the ‘brms’ package [version 2.22.0; (63)], which fits multilevel Bayesian models using the probabilistic programming language, ‘Stan’ (71), accessed via the ‘rstan’ package (version 2.32.7 (68, 69)). All models employed 4 chains with 8,000 iterations (including 2,000 and 6,000 warm-up and sampling iterations, respectively). Diagnostic checks included checks of MCMC performance and also several prior and posterior validation checks (Supplementary File 3), with full details of all statistical analyses (including the code used and all statistical output) available online: https://github.com/AliG71/hair_cortisol. To verify MCMC performance, models were checked for convergence, by assessing R-hat values and inspecting both trace and trace-rank plots (Supplementary File 3), whilst resolution was assessed by calculating effective sample sizes (Table 4). In all final models, effective sample sizes were always acceptable (typically 10,000–20,000 depending upon the model).

Table 4

Table 4. Summary of causal associations for the final models where log hair cortisol concentration was the outcome variable.

Verifications of the suitability of prior probabilities included initial graphical modelling, to simulate the expected shape of the distribution, prior predictive simulations and power scaling sensitivity analysis [powerscale_sensitivity function of the ‘priorsense’ package (67)]. Verification checks on the posterior distributions included a visual inspection of a pairs plot, a graphical posterior predictive check, graphical comparisons of individual draws against the observed data for each predictor variable and several other graphical checks (Supplementary File 3).

Finally, before testing on the actual study data, models were checked using simulated data to ensure that model fitting had worked correctly. For this, different simulated datasets of 200 visits [100 each for the first (V0) and second (V1) visits from 100 dogs] were created, each specifying different effect sizes for the causal predictors of each model (e.g., in the body fat model, it was assumed that the beta coefficient for the causal effect could be of 0.0, 0.2 or 0.5). All model estimates for the beta coefficient reliably reproduced what was expected in the simulated dataset (https://github.com/AliG71/hair_cortisol).

Analyses of the posterior distribution

Posterior density distributions for all parameters in the final models were calculated as described above and summarised using means (for central tendency) and 97% highest posterior density intervals (97% HPDI, for limits of the credible interval), unless otherwise indicated. Graphical visualisation of posterior probability distributions included plots of posterior densities (displaying highest density intervals), plots of conditional effects and plots comparing the prior and posterior probability distributions (to illustrate the relative contribution of the prior and sample data). Hypothesis tests were conducted to determine the probability that the effects of each causal predictor were positive (or negative), and posterior predictions were used (from new simulated data) to estimate causal associations. Finally, the Bayes R² metric was calculated using the bayes_R2 function of the ‘rstantools’ package (69), which is similar to a conventional R² statistic from least-squares linear regression, with the results “quantifying the fit of the model to the data at hand” (69).

Sensitivity analyses

Sensitivity analyses included testing the sensitivity of the posterior distribution to the choice of prior distribution, as described above. Other analyses were conducted, on the body fat and comorbidity models, to determine the sensitivity of posterior distribution estimates to measurement error (given known variability in the hair cortisol assay), season of sampling, missing data (body fat model only) and setting a neutral regularising prior for comorbidity (comorbidity model only). Measurement error in the response variable (logHCC) was accommodated using the mi() syntax in the ‘brms’ package (63), and assuming that the coefficient of variability of hair cortisol measurement was 11.8% (37). Although not on the causal pathway for either model (Figure 1; Supplementary File 1), a sensitivity analysis was conducted to determine the effect of sample season on causal estimates in the body fat and comorbidity models. This was undertaken by recomputing posterior probability distribution after adding the sampling season variable to each model. To assess a possible effect of missing data in the body fat variable, the mi() syntax from the ‘brms’ package (63) was again used; this involved fitting a multivariate Bayesian multilevel model, whereby the missing data for body fat and log hair cortisol are simultaneously predicted (see the sensitivity analyses report in the online material at: https://github.com/AliG71/hair_cortisol). Using this approach, missing data are handled as additional parameters for estimation from a single joint posterior distribution. Knowledge of variables causally-associated with the missing variable can also be used to inform the imputation process. Together, this leads to more accurate and honest credible intervals than other imputation approaches such as single imputation. For the body fat model, a Student’s t likelihood function was chosen, and incorporated appropriate predictor variables as indicated by the DAG (Figure 1; Supplementary File 1; e.g., sex, age, breed group and comorbidity). Finally, to determine whether the final choice of prior for comorbidity (marginally-positive, weakly-regularising; mean 0.25, sigma 1) had unduly influenced the posterior probability density, the comorbidity model was rerun with a different prior that was neutral and weakly regularising (mean 0, sigma 1).

Ethics and welfare considerations

The study has been reported in accordance with the Animal Research: Reporting of In Vivo Experiments (ARRIVE) guidelines.¹ The study received approval from both the University of Liverpool Veterinary Research Ethics Committee (RETH000353 and VREC793) and the Royal Canin Ethical Review Committee (150720–55). All owners gave informed, written consent allowing their dog to participate. Clinical procedures complied with relevant guidelines (e.g., standard operating procedures) and regulations. Foods used were commercially available therapeutic diets commonly used by veterinarians to manage obesity, and given for the clinical benefit of the study dogs. Neither the clinical procedures used nor the clinical use of the therapeutic diets were deemed to involve animal experimentation, falling outside the remit of national legislation (e.g., the revised Animals [Scientific Procedures] Act 1986).

Results

Characteristics of dogs, hair samples and HCC results

In total, 73 hair samples were obtained, comprising paired samples from 21 dogs (before [V0] and after [V1] weight reduction), and single samples (before weight reduction only [V0]) from a further 31 dogs. Full details of all baseline characteristics are shown in Table 3. Forty-four of the 52 dogs (85%) had one or more comorbidities (median 1, range 0–3). Twenty-five dogs ate the dry therapeutic food during their period of weight reduction (16 dogs with paired samples; 9 dogs with single samples), whilst the remaining 27 consumed a combination of wet and dry therapeutic food (12 dogs with paired samples; 15 dogs with single samples). In the dogs with paired samples, mean weight loss was 26% (SD 6.8%) of starting weight, at a rate of 0.7% (SD 0.33) per week, with median body fat mass being 42% (range 31–61%) and 33% (range 18–50%) before and after weight reduction, respectively.

Table 3

Table 3. Baseline variables in dogs in the study.

There were various reasons why 31 of the dogs only contributed a single hair sample. Just over half (16 dogs) had reasons related to the COVID-19 pandemic: of these, 3 dogs reached target weight during the pandemic, but a face-to-face follow-up visit was not possible; the remaining 13 pandemic-affected dogs were lost to follow-up, mostly because owners had found it difficult to implement a therapeutic weight reduction protocol. Of the 15 dogs not affected by the COVID-19 pandemic, 3 dogs completed therapeutic reduction but their owners decided not to return for the follow up; 2 dogs were euthanased for other reasons (e.g., metastatic pulmonary adenocarcinoma, old age) before reaching target weight; and the remaining 10 dogs were lost to follow-up (4 due to poor compliance; 4 stopped responding to communications; 1 moved away from the area; 1 had transport issues).

Hair cortisol concentrations

Details of HCC results are shown in Table 3. The mean HCC for all samples was 10.4 (SD 19.52) pg/mg (logHCC 1.3, SD 1.36). In the dogs providing single samples (before therapeutic weight reduction, V0), mean HCC was 11.5 (SD 21.91) pg/mg (log HCC 1.4, SD 1.33) whereas, in those providing paired samples, HCC was 7.6 (SD 14.29) pg/mg (log HCC 1.1, SD 1.23) before (V0) and 3.6 (SD 3.74) pg/mg (log HCC 0.9 SD 0.94) after (V1) weight reduction.

Causal associations between hair cortisol and different explanatory variables

The results of estimates of the causal effects from each model are summarised in Table 4 and Supplementary File 5, with full details of each statistical model and validation available online: https://github.com/AliG71/hair_cortisol.

Age model

Based on the DAG (Figure 1; Supplementary File 1), the only adjustment variable required for the age model was breed group. Although the probability distribution for the effect of age (beta coefficient) spanned zero (Supplementary File 5), and the majority was negative (mean −0.12 per SD; 97% HPDI -0.33, 0.10). A weakly-negative trend was seen across the age range (Supplementary File 5), equating to an average decrease in HCC of ~1.6 pg/mg (~0.1 logHCC) for each 1-unit (~36 months) increase in standardised age. Therefore, we estimated an 88% probability of there being a negative effect of age on logHCC.

Sex model

Based on the DAG (Figure 1; Supplementary File 1), no adjustment variables were required for the sex model. The probability distribution for the effect of sex again spanned zero (Supplementary File 5), but did not clearly favour either a positive or negative predictions (mean 0.14, 97% HPDI -0.29, 0.60). Therefore, we estimated a 75% probability of there being a positive causal association between sex and logHCC.

Breed group model

Based on the DAG (Figure 1; Supplementary File 1), no adjustment variables were again required for the breed group model. Compared with mixed breed as the reference category, probability distributions for all other categories spanned zero (Supplementary File 5), and did not suggest any clear positive or negative associations with logHCC (CKCS mean 0.14, 97% HPDI -0.72, 0.92; pug mean 0.19, 97% HPDI -0.59, 0.62; retriever mean 0.09, 97% HPDI -0.52, 0.70; other breed mean 0.10, 97% HPDI -0.46, 0.65). Therefore, the probability of positive causal associations with logHCC were estimated to be 66%, 71%, 63% and 65% for the CKCS, pug, retriever and other breed groups, respectively.

Coat colour model

Based on the DAG (Figure 1; Supplementary File 1), the only adjustment variable required for the coat colour model was breed group. Compared with light coat colour as the reference category, the probability distributions for mixed and dark coat colour spanned zero (Supplementary File 5) but were predominantly negative in both instances (mixed: mean −0.50, 97% HPDI -1.13, 0.13; dark: mean −0.33, 97% HPDI -0.83, 0.17). Therefore, we estimated 93% probability and 96% probability that logHCC is less in dogs with dark or mixed coat colour, respectively than in dogs with light coat colour.

Season of sampling model

Based on the DAG (Figure 1; Supplementary File 1), no adjustment variables were again required for the sampling season model. Compared with spring as the reference category, all probability distributions spanned zero (Supplementary File 5), although those for summer (mean −0.18; 97% HPDI -0.73, 0.37) and autumn (mean −0.32; 97% HPDI -0.89, 0.29) were majority negative, whilst that for winter was majority positive (mean 0.26; 97% HPDI -0.34, 0.88). Therefore, the probability of negative causal associations between summer and autumn and logHCC were 77% and 88%, respectively, whilst the probability of a positive causal association between winter and log HCC was 83%.

Comorbidity model

Based on the DAG (Figure 1; Supplementary File 1), the adjustment set required for the comorbidity model included age, sex and breed group. The entire posterior probability distribution was positive (mean 0.65; 97% HPDI: 0.11, 1.16; Bayes R² 0.13; Figure 5a), meaning that having a comorbidity was associated with an average increase in HCC of 10.4 pg/mg, compared with no comorbidity (Figure 5b). The relative contributions of the prior expectations and the observed data were determined in two ways; first, the prior and posterior density distributions of the beta coefficient for comorbidity were visually compared (Figure 6a); second, a visual comparison was made of 50 random draws taken from the prior (Figure 6b) and posterior (Figure 6c) probability distributions. Prior to observing the data, the model expected the causal association between comorbidity and logHCC to be slightly positive, on average, but with wide degree of uncertainty as indicated by the mean effect of the prior being >0 but with wide variability (Figure 6a) and the fact that individual posterior predictions of logHCC were slightly more likely to be positive than negative, albeit to varying degrees (Figure 6b). After seeing the data, the model was both confident that the causal association would be positive, and relatively confident about its magnitude, as indicated by the narrow range of the posterior distribution (Figure 6a) and the close clustering of the posterior predictions of logHCC (Figure 6c). Therefore, we estimated a > 99% probability of a positive causal association between comorbidity and logHCC.

Figure 5

Figure 5. (a) Highest probability density interval (a.k.a. highest density interval) plot for the causal association between comorbidity and log hair cortisol concentration (logHCC). The y-axis depicts the density distribution, whilst the x-axis depicts possible values for the beta parameter of the causal association. Different probability intervals are depicted by colour (red 65%; purple 70%; orange 80%; yellow 89%; green 97%; blue 100%). The feint vertical dotted line intersects the x-axis at zero, indicating the point where the beta coefficient would be neither positive of negative. Most of the posterior probability density is greater than zero, indicating a > 99% probability for a positive causal association. (b) Conditional effect for the association between comorbidity and logHCC, generated using the posterior_predict function of the ‘brms’ package [version 2.22.0; (62)]. The points represent estimates, whilst the lines represent the 97% compatibility interval (97%-CI). This method estimates uncertainty in predictions from both the statistical model and residual error, which better reflects the complete range of plausible outcomes of the scientific model. Mean predicted logHCC is greater in dogs with comorbidities, with a wider uncertainty range, than in dogs without comorbidities.

Figure 6

Figure 6. A comparison of prior and posterior distributions for the causal association between comorbidity and log hair cortisol concentrations (logHCC). (a) Comparison of prior and posterior probability densities (thin blue lines). Thick vertical lines within these densities represent the mean and the shaded blue regions represent the central (50%) compatibility interval. Differences in in the shape and position of the posterior distributions, relative to the prior distribution, indicate the effect of the observed data on the predictions. (b,c) Visual comparison of predicted differences in logHCC between dogs with and without comorbidities, based on individual predictions created using 50 random draws taken from the prior (b) or posterior (c) probability distribution. Data are displayed as combined box-and-whisker and point plots, with thick horizontal lines representing the median, and upper and lower limits of the boxes represent the inter-quartile range (IQR). The upper and lower whiskers extend as far as the largest or smallest, respectively, values that are no further than 1.5 × IQR from the IQR, whilst outlying points are shown as small black dots. Prior to observing the data, the model predicts a slightly positive average causal association, but with wide degree of uncertainty, as indicated by the wide density distribution spanning zero (a), and only a marginal average difference between the box and whisker plots for each group (b). After seeing the data, the model is confident that the causal association between comorbidity and logHCC is positive, as indicated by the narrower density distribution almost-completely >0, and greater difference between box-and-whisker plots (C).

Sensitivity analyses (https://github.com/AliG71/hair_cortisol) demonstrated that these results were robust to the effects of measurement error (mean 0.65, 97% HPDI 0.11, 1.17; probability >99%; 10.4 pg/mg mean increase in HCC when one or more comorbidity present), sampling season (mean 0.60, 97% HPDI 0.03, 1.16; probability >99; 9.6 pg/mg mean increase in HCC when one or more comorbidity present) and using a neutral, weakly-regularising prior for comorbidity (mean 0.63, 97% HPDI 0.10, 1.14; probability >99%; 10.1 pg/mg mean increase in HCC when one or more comorbidity present).

Body fat model

Based on the DAG (Figure 1; Supplementary File 1), the adjustment set required to estimate the causal association between body fat and logHCC, included age, sex, breed group and comorbidity. The vast majority of the probability distribution for the beta coefficient was positive (mean 0.33 per SD; 97% HPDI 0.03, 0.59; Bayes R² 0.18; Figure 7a), and a positive trend was seen across the body fat range (Figure 7b), equating to an average increase in HCC of 5.3 pg/mg for each increase of 1 unit in standardised fat (~10% body fat mass). The relative contributions of our prior expectations and the observed data were assessed in two ways; first, the prior and posterior density distributions of the beta coefficient for body fat from the final model were visually compared (Figure 8a); second, a visual comparison was made of the individual regression lines predicted from 50 random draws taken from the prior (Figure 8b) and posterior (Figure 8c) probability distributions. Prior to observing the data, the model was both neutral and uncertain about the possible effect of body fat on logHCC, as indicated by the mean effect of the prior being 0 (Figure 8a) and the fact that individual regression lines could be either positive or negative to varying degrees (Figure 8b). After seeing the data, the model was both confident of a positive causal association, and relatively confident about its magnitude, as indicated by the narrow range of the posterior distribution (Figure 8a) and the close clustering of the regression line slopes (Figure 8c). Therefore, we estimated a 99% probability of a positive causal association between body fat and logHCC. The conditional effect of body fat on logHCC, stratified by comorbidity and compared with the observed data is shown in Figure 9.

Figure 7

Figure 7. (a) Highest probability density interval (a.k.a. highest density interval) plot for the causal association between body fat and log hair cortisol concentration (logHCC). The y-axis depicts the density distribution, whilst the x-axis depicts possible values for the beta parameter of the causal effect. Different probability intervals are depicted by colour (red 65%; purple 70%; orange 80%; yellow 89%; green 97%; blue 100%). The feint vertical dotted line intersects the x-axis at zero, indicating the point where the beta coefficient would be neither positive of negative. Most of the posterior probability density is greater than zero, indicating a 98% probability that the causal association is positive. (b) Conditional effect for the causal association between body fat and logHCC. The blue line represents the mean estimate for logHCC, across the body fat range, whilst the shaded region represents the 95% credible interval. Estimates were generated using the posterior_predict function of the ‘brms’ package [version 2.22.0 (62)], which returns the posterior mean and 95% credible interval for each data point, thereby incorporating uncertainty in predictions from the posterior distribution (model predictions) and uncertainty due to residual error (from individual data points). Credible intervals generated in this way reflect the complete range of plausible outcomes from the scientific model. Overall, there is a positive linear relationship between body fat and logHCC.

Figure 8

Figure 8. A comparison of prior and posterior distributions for the causal association between body fat and log hair cortisol concentrations (logHCC). (a) Comparison of prior and posterior probability densities (thin blue lines). Thick vertical lines within these densities represent the mean and the shaded blue regions represent the central (50%) compatibility interval. Differences in the shape and position of the posterior distributions, relative to the prior distribution, indicate the effect of the observed data on the predictions. (b,c) Visual comparison of predicted differences in regression lines. Each line is an individual prediction of the regression slope, calculated from the intercept and beta coefficient for body fat, using one of 50 random draws from the prior (b) or posterior (c) probability distributions. Prior to observing the data, the model is both neutral and uncertain about the possible causal association between body fat and logHCC, as indicated by the wide density distribution spanning zero (a) and wide range of possible regression slopes (b), albeit limited to a physiologically plausible range by the regularising nature of the prior. After seeing the data, the model is more that the association is positive, as indicated by a narrower density distribution that is predominantly above zero, and a more limited range of regression slopes, all of which are positive.

Figure 9

Figure 9. Conditional effect of body fat on log hair cortisol concentration (logHCC), stratified by comorbidity (red: no comorbidity, blue: comorbidity) and compared with the observed data. Solid lines represent mean predictions in logHCC across the body fat range, with shaded areas representing 95% compatibility intervals, and points representing observed results from individual dogs. Estimates were generated using the posterior_epred function of the ‘brms’ package [version 2.22.0 (62)], which only incorporates uncertainty from the posterior distribution (model predictions) and, therefore, uncertainty bands are narrower than those seen in Figure 7C. An equivalent plot using posterior_predict is included in the Supplementary File 5. Again, a positive causal association between body fat and logHCC is evident, with values being greater and more variable in dogs with comorbidities.

Sensitivity analyses (Supplementary File 5) demonstrated that these results were relatively robust to the effects of measurement error (mean 0.32, 97% HPDI 0.03, 0.60; probability 98%; 5.1 pg/mg mean increase in HCC per 10% body fat mass increase), season of sampling (mean 0.28, 97% HPDI -0.11, 0.62; probability 94%; 4.5 pg/mg mean increase in HCC per 10% body fat mass increase) and missing data (mean 0.27, 97% HPDI 0.01, 0.55; probability 98%; 4.3 pg/mg mean increase in HCC per 10% body fat mass increase.

Reverse causality model

Although we had assumed that changes in body fat mass would cause changes in HPA function, the possibility of reverse causality, namely that an upregulated HPA would cause increased body fat mass, could not be discounted. To explore this possibility, the DAG was redrawn with HCC and body fat as exposure and outcome, respectively (Supplementary File 1). The same adjustment set as with the body fat model (e.g., age, sex, breed group and comorbidity) was required for this causal analysis, which was tested twice, first using all study data (model 1, Table 5; Figures 10a,c) and also using data from visit 0 only (model 2, Table 5; Figures 10b,d). In both models, posterior probability densities were wide, spanning zero (Figures 10a,b), albeit with slightly more of the probability mass being positive (all data 70%; visit 0 data only 87%). Further, the linear trend was relatively flat across the logHCC range (Figures 10c,d), with the credible interval being broad and including both positive and negative slopes. Such results suggest uncertainty in the estimate of this effect, and would not be consistent with a there being a convincing reverse causality effect.

Figure 10

Figure 10. (a,b) Highest probability density intervals (a.k.a. highest density interval) plot for the causal association between log hair cortisol concentration (logHCC) and body fat, taken from the reverse causality models and either using all study data (a) or only data from visit 0 (b). The y-axis depicts the density distribution, whilst the x-axis depicts possible values for the beta parameter of the causal effect. Different probability intervals are depicted by colour (red 65%; purple 70%; orange 80%; yellow 89%; green 97%; blue 100%). The feint vertical dotted line intersects the x-axis at zero, indicating the point where the beta coefficient would be neither positive of negative. In both models, the posterior probability density spans zero, indicating that the effect could credibly be either positive or negative, albeit with positive effects being more likely given that 70% (a: all data) or 87% (b: visit 0 data) of the probability density is positive. (c,d) Conditional effects for the causal association between logHCC and body fat taken from the reverse causality models and either using all study data (c) or only data from visit 0 (d). The blue line represents the mean estimate for logHCC, across the body fat range, whilst the shaded region represents the 95% credible interval. Estimates were generated using the posterior_predict function of the ‘brms’ package [version 2.22.0; (62)], which returns the posterior mean and 95% credible interval for each data point, thereby incorporating uncertainty in predictions from the posterior distribution (model predictions) and uncertainty due to residual error (from individual data points). Credible intervals generated in this way reflect the complete range of plausible outcomes from the scientific model. Overall, the linear relationship between logHCC and body fat is relatively flat, with a broad credible interval that includes horizontal).

Table 5

Table 5. Summary of causal associations for the final models where log hair cortisol concentration was the outcome variable.

Discussion

The objective of the present study was to use a Bayesian workflow to estimate possible causal associations between HCC and several different variables in dogs with obesity. The variables for which causal associations were most plausible were body fat and having at least one comorbidity, given the overwhelmingly positive posterior probability densities obtained in the final models. The reasons for a possible positive causal association between body fat and HCC are not clear but a similar association is seen in humans (72). Of course, the causal association between obesity and HPA upregulation might flow in the opposite direction, so-called reverse causation, where upregulated HPA causes an increased body fat mass. Evidence for this includes the fact that cortisol promotes adipose tissue redistribution (to the abdominal region), increases appetite and, in humans at least, promotes a preference for foods of greater palatability (73), all of which could lead to adipose tissue gain. Increased cortisol concentration is also observed in people who gain weight because of stress (74). To address this, we created a final causal model with body fat and logHCC as outcome and causal predictors, respectively, and utilising the same adjustment variables as required for the body fat model (Supplementary File 1). We chose to assess this association both using all study data and only data from visit 0. We included the second analysis because we were concerned that therapeutic weight reduction might affect HCC, as reported in a previous study in dogs (20), thereby obscuring any causal effect of HPA upregulation on body fat mass. In contrast to the result of the body fat model, where the probability distribution for a causal effect was overwhelmingly positive (99%), the probability density spanned zero in both the reverse causality models, i.e., not consistent with a causal effect in this direction. Taken together, these results suggest that, if a causal effect exists, it is most likely to be in the direction of body fat causing HPA upregulation.

For dogs where paired samples were available, mean HCC was 7.6 pg/mg and 3.6 pg/mg in samples taken before and after therapeutic weight reduction, respectively. A major study limitation was the fact that it was not possible to create an unbiased statistical model for the effect of therapeutic weight reduction because, based on our DAG, we could not correct for possible confounding from unmeasured variables. Such unmeasured variables could include owner factors (e.g., owner attitudes and behaviours in implementing a therapeutic weight reduction protocol), environmental factors (e.g., the living environment of the dog) and, of course, the impact of the COVID-19 pandemic; all such variables could plausibly affect both the body fat percentage of dogs and the confounding variables (e.g., affecting attendance of owners at visits and whether an individual dog would successfully reach its target). Considering the findings of the current study, a prospective study should now be considered where the effect of therapeutic weight reduction on HCC could be assessed prospectively.

As well as a possible direct effect of excessive adipose tissue mass on HPA upregulation, obesity might have an indirect effect either by predisposing to or exacerbating other diseases; indeed, many possible disease associations have been identified in previous research (11). In the current study, most dogs had at least one comorbidity, with some having multiple comorbidities. The estimated causal association between comorbidity and HCC was highly variable; not only was HCC greater by an average of 10.4 pg/mg in dogs with comorbidities, but concentrations within the group were much more variable than in the dogs without comorbidities. These findings are, perhaps, not surprising given the many different comorbidities present, with differing severity, and the fact that the pathogenetic mechanisms are likely to be different amongst different diseases. For example, orthopaedic diseases were particularly common in the current study; such diseases might activate stress pathways because of chronic pain, trauma to local tissues caused by the injury or subsequent surgery, inflammation, and by locally triggering endogenous corticosteroid production, as seen in human studies (75–77). Further studies would be required to confirm the mechanisms by which excess adipose tissue leads to increased HCC in dogs.

Possible causal associations were also suggested for some of the other variables, including age and hair coat although posterior probability distributions indicated that such associations were less certain. In this respect, the probability of there being a negative association between age and HCC was estimated to be 88%, whilst the probability of there being negative associations between mixed or dark coat colour and HCC were estimated to be 96 and 93%, respectively. A possible association with coat colour was suggested in previous research (8), and this was the justification for setting a marginally-negative regularising prior in this model. The reason for this association is not fully understood but is likely due to hair-related rather than systemic factors because differences in salivary cortisol are not seen in dogs with different coat colour (8). This has resulted in suggestions that hair of different colour might sequester cortisol to differing degrees, perhaps associated with differences in the size of melanin granules (78). In contrast to the association with hair colour, age has not previously been identified as being associated with HCC in dogs (8), whilst the effect of age in humans is non-linear, being greater in children and elderly, and lowest in middle age (79). Only adult dogs were assessed in the current study (with an age range at the initial visit of 2 to almost 14 years). Therefore, further research would be required, including sampling from growing dogs, to explore more fully a possible causal association between age and HCC in dogs.

There are several limitations that should be considered:

1. The study population was relatively small, with considerable variability in baseline attributes amongst dogs including the presence of comorbidities, and we did not control for the therapy that dogs received. This will have created noise within the dataset, which might have masked the actual causal associations with the variables studied.

2. Having missing data in the body fat percentage variable created challenges for data analysis, although our sensitivity analyses suggested that the impact of this was minimal.

3. It was not possible to examine the causal association between weight loss success and HCC, for example, by creating a binary variable classifying dogs into those that completed and those that ended their therapeutic weight reduction prematurely. Based on our scientific model (Figure 1; Supplementary File 1), it was not possible to identify a set of adjustment variables that allowed this causal association to be estimated. Rather than reporting an erroneous result, with potentially misleading conclusions, we preferred not to attempt such an analysis.

4. It was also not possible to develop a causal model for an effect of visit (before vs. after weight reduction; Figure 1; Supplementary File 1). Therefore, questions of the effect of visit and non-compliance with weight reduction are arguably best answered using a prospective study with an appropriate design.

5. All statistical analyses were based on a scientific model, codified as a DAG, and there might have been errors and omissions, which adversely impact the accuracy of any causal effects.

6. A limitation of the DAG method is that causal pathways typically need to be one-way and cyclical paths are not allowed (i.e., one variable cannot influence other variables which, in turn, influence the first variable). Possible inverse causality in the relationship between changes in body fat and HPA function was considered, as discussed above. However, the direction of association between body fat and having comorbidities is a similar consideration. In our scientific model, we assumed that comorbidities would have an impact on body fat percentage, but the association might be reversed, at least for some diseases. We chose the former rather than the latter on the basis that data on comorbidities were cross-sectional, in that all comorbidities were present at the time of the initial visit. Therefore, the effects of any such comorbidity on body fat mass, for example by affecting food intake or physical activity, would already be evident. In contrast, to test the causal effect of body fat mass on comorbidities, a longer-term cohort study would be needed whereby the effect of initial body fat mass on the future development of a comorbidity could be assessed, as with a recent study examining the association between obesity and future diabetes mellitus (80). Since we did not have body fat data prior to the development of the comorbidities, it would be difficult to study such a causal effect. Further, any monitoring period in such a study period would need to be of sufficient duration to maximise the chances of fully capturing the effect; indeed, a monitoring period of at least 4 years was assessed in the previous study of obesity and DM (80).

Besides these limitations, the findings of the current study should also be interpreted considering possible methodological limitations inherent to hair cortisol analysis. First, HCC can vary with body region and cortisol deposition may not be uniform; in a study involving both humans and other animals, HCC could differ by >20% in hair samples taken from different locations (e.g., head vs. limbs) (81). We controlled for this by always sampling from the same two regions and pooling the harvested hair. These sites were chosen for reasons of convenience since they happened to be regions where clipping was already required for medical procedures (e.g., blood sampling and intravenous catheterisation). Although this approach was favoured for ethical and welfare reasons, our results might not be directly comparable to other studies where different body regions were sampled.

We also did not account for possible variability caused by differences in hair growth rates or the stage of the hair growth cycle, not least since dog hair growth is often cyclical and can be influenced by season (5). Compared with spring, the probability of negative causal associations between summer or autumn and HCC were 77 and 88%, respectively, whilst the probability of a positive causal association between winter and HCC was 83%. Given such probabilities, seasonal effects are certainly possible, but by no means probable. In our scientific model (Figure 1), back door pathways between season of sampling and either body fat percentage or comorbidly were thought to be implausible and, therefore, sampling season was not likely to be a confounder for either variable. These assumptions were confirmed in sensitivity analyses whereby the effects of both body fat percentage and comorbidity were largely unchanged when sampling season was included in modelling.

Finally, there is some evidence for local cortisol production in the skin and hair; for example, in guinea pigs, systemically-administered, radiolabelled cortisol accounted for only a small fraction of the cortisol found in hair, with the remainder probably arising from local follicular synthesis (82). Such local production of cortisol or cortisol-like compounds might also contribute to HCC in dogs, as suggested for dogs with HAC (5). Therefore, factors affecting the skin (e.g., skin inflammation, hyperpigmentation, or topical steroid exposure) might feasibly affect HCC independently of systemic cortisol status.

In conclusion, increased body fat and the presence of one or more comorbidities are causally associated with increased HCC in dogs; prospective studies should assess the impact of therapeutic weight reduction.

Data availability statement

All data generated or analysed during this study are included in this published article and its supplementary information files (Supplementary File 2). The data, statistical code and all statistical reports are also available online at: https://github.com/AliG71/hair_cortisol.

Ethics statement

The animal studies were approved by University of Liverpool Veterinary Research Ethics Committee and The Royal Canin Ethical Review Committee. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent was obtained from the owners for the participation of their animals in this study.

Author contributions

KT: Data curation, Investigation, Visualization, Writing – original draft, Writing – review & editing. GW-L: Data curation, Methodology, Project administration, Writing – review & editing. JF: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Validation, Writing – review & editing. XL: Validation, Writing – review & editing. AG: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. The study was funded by Royal Canin.

Acknowledgments

The authors wish to acknowledge the referring veterinarians for referring cases and the owners of all dogs for allowing them to participate. The authors also thank Vincent Biourge for scientific comments. Preliminary results were presented as an abstract at the European College of Veterinary Internal Medicine Congress, Maastricht, Netherlands, Sept 2025.

Conflict of interest

The study was funded by a grant from Royal Canin, a division of Mars Petcare. The funder had the following involvement in the study: the company manufactured the diets fed in this study, provided support to the infrastructure of the Liverpool obesity care clinic, which includes equipment funding, clinical diagnostics, clinical consumables, and therapeutic food. JF and XL are employees of Royal Canin. AG and GW-L are employees of the University of Liverpool but their positions are funded by Royal Canin. Both have received financial remuneration and gifts for providing educational material, speaking at conferences, and consultancy work.

The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The authors declare that no Gen AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2025.1695345/full#supplementary-material

Footnotes

1. ^https://arriveguidelines.org/arrive-guidelines

References

1. Meyer, JS, and Novak, MA. Minireview: hair cortisol: a novel biomarker of hypothalamic-pituitary-adrenocortical activity. Endocrinology. (2012) 153:4120–7. doi: 10.1210/en.2012-1226

PubMed Abstract | Crossref Full Text | Google Scholar

2. Corradini, S, Accorsi, PA, Boari, A, Beghelli, V, Mattioli, M, Famigli-Bergamini, P, et al. Evaluation of hair cortisol in the diagnosis of hypercortisolism in dogs. J Vet Intern Med. (2013) 27:1268–72. doi: 10.1111/jvim.12135

PubMed Abstract | Crossref Full Text | Google Scholar

3. Mack, Z, and Fokidis, HB. A novel method for assessing chronic cortisol concentrations in dogs using the nail as a source. Domest Anim Endocrinol. (2017) 59:53–7. doi: 10.1016/j.domaniend.2016.11.003

PubMed Abstract | Crossref Full Text | Google Scholar

4. Roth, LS, Faresjö, Å, Theodorsson, E, and Jensen, P. Hair cortisol varies with season and lifestyle and relates to human interactions in German shepherd dogs. Sci Rep. (2016) 6:19631. doi: 10.1038/srep19631

PubMed Abstract | Crossref Full Text | Google Scholar

5. Ouschan, C, Kuchar, A, and Möstl, E. Measurement of cortisol in dog hair: a noninvasive tool for the diagnosis of hypercortisolism. Vet Dermatol. (2013) 24:e93-4:428–31. doi: 10.1111/vde.12043

Crossref Full Text | Google Scholar

6. Hodes, A, Meyer, J, Lodish, MB, Stratakis, CA, and Zilbermint, M. Mini-review of hair cortisol concentration for evaluation of Cushing syndrome. Expert Rev Endocrinol Metab. (2018) 13:225–31. doi: 10.1080/17446651.2018.1517043

PubMed Abstract | Crossref Full Text | Google Scholar

7. Grigg, EK, Nibblett, BM, Robinson, JQ, and Smits, JE. Evaluating pair versus solitary housing in kennelled domestic dogs (Canis familiaris) using behaviour and hair cortisol: a pilot study. Vet Rec Open. (2017) 4:e000193. doi: 10.1136/vetreco-2016-000193

PubMed Abstract | Crossref Full Text | Google Scholar

8. Bowland, GB, Bernstein, RM, Koster, J, Fiorello, C, Brenn-White, M, Liu, J, et al. Fur color and nutritional status predict hair cortisol concentrations of dogs in Nicaragua. Front Vet Sci. (2020) 7:565346. doi: 10.3389/fvets.2020.565346

PubMed Abstract | Crossref Full Text | Google Scholar

9. Knezevic, E, Nenic, K, Milanovic, V, and Knezevic, NN. The role of cortisol in chronic stress, neurodegenerative diseases, and psychological disorders. Cells. (2023) 12:2726. doi: 10.3390/cells12232726

PubMed Abstract | Crossref Full Text | Google Scholar

10. Alotiby, A. Immunology of stress: a review article. J Clin Med. (2024) 13:6394. doi: 10.3390/jcm13216394

PubMed Abstract | Crossref Full Text | Google Scholar

11. Lund, EM, Armstrong, J, Kirk, CA, and Klausner, JF. Prevalence and risk factors for obesity in adult dogs from private us veterinary practices. Int J Appl Res Vet Med. (2006) 4:177–86. https://jarvm.com/articles/Vol4Iss2/Lund.pdf

Google Scholar

12. German, AJ, Holden, SL, Wiseman-Orr, ML, Reid, J, Nolan, AM, Biourge, V, et al. Quality of life is reduced in obese dogs but improves after successful weight loss. Vet J. (2012) 192:428–34. doi: 10.1016/j.tvjl.2011.09.015

PubMed Abstract | Crossref Full Text | Google Scholar

13. Salt, C, Morris, PJ, Wilson, D, Lund, EM, and German, AJ. Association between life span and body condition in neutered client-owned dogs. J Vet Intern Med. (2019) 33:89–99. doi: 10.1111/jvim.15367

PubMed Abstract | Crossref Full Text | Google Scholar

14. Henegar, JR, Bigler, SA, Henegar, LK, Tyagi, SC, and Hall, JE. Functional and structural changes in the kidney in the early stages of obesity. J Am Soc Nephrol. (2001) 12:1211–7. doi: 10.1681/ASN.V1261211

PubMed Abstract | Crossref Full Text | Google Scholar

15. Marshall, WG, Hazewinkel, HA, Mullen, D, De Meyer, G, Baert, K, and Carmichael, S. The effect of weight loss on lameness in obese dogs with osteoarthritis. Vet Res Commun. (2010) 34:241–53. doi: 10.1007/s11259-010-9348-7

PubMed Abstract | Crossref Full Text | Google Scholar

16. Tvarijonaviciute, A, Ceron, JJ, Holden, SL, Cuthbertson, DJ, Biourge, V, Morris, PJ, et al. Obesity-related metabolic dysfunction in dogs: a comparison with human metabolic syndrome. BMC Vet Res. (2012) 8:147. doi: 10.1186/1746-6148-8-147

PubMed Abstract | Crossref Full Text | Google Scholar

17. Tvarijonaviciute, A, Ceron, JJ, Holden, SL, Biourge, V, Morris, PJ, and German, AJ. Effect of weight loss in obese dogs on indicators of renal function or disease. J Vet Intern Med. (2013) 27:31–8. doi: 10.1111/jvim.12029

PubMed Abstract | Crossref Full Text | Google Scholar

18. Wardle, J, Chida, Y, Gibson, EL, Whitaker, KL, and Steptoe, A. Stress and adiposity: a meta-analysis of longitudinal studies. Obesity (Silver Spring). (2011) 19:771–8. doi: 10.1038/oby.2010.241

PubMed Abstract | Crossref Full Text | Google Scholar

19. Larsen, SC, Turicchi, J, Christensen, GL, Larsen, CS, Jørgensen, NR, Mikkelsen, MK, et al. Hair cortisol concentration, weight loss maintenance and body weight variability: a prospective study based on data from the European NoHoW trial. Front Endocrinol (Lausanne). (2021) 12:655197. doi: 10.3389/fendo.2021.655197

Crossref Full Text | Google Scholar

20. Kim, K, Song, B, Kim, D, Kim, DH, Lee, HJ, and Kim, G. Effect of leash walking on weight loss and assessment of hair cortisol in overweight dogs. Comp Exerc Physiol. (2024) 20:283–91. doi: 10.1163/17552559-20231019

Crossref Full Text | Google Scholar

21. Pearl, J. Causal inference in statistics: an overview. Stat Surv. (2009) 3:96–146. doi: 10.1214/09-SS057

Crossref Full Text | Google Scholar

22. Donati, PA, Tunesi, M, Portela, DA, Maxwell, EA, Campana, JP, Reina, JP, et al. Preoperative septic peritonitis, hypotension, and reason for surgery are risk factors for small intestine dehiscence in dogs: a directed acyclic graph approach. J Am Vet Med Assoc. (2025) 263:1–8. doi: 10.2460/javma.24.12.0791

PubMed Abstract | Crossref Full Text | Google Scholar

23. Weng, HY, Kass, PH, Hart, LA, and Chomel, BB. Risk factors for unsuccessful dog ownership: an epidemiologic study in Taiwan. Prev Vet Med. (2006) 77:82–95. doi: 10.1016/j.prevetmed.2006.06.004

PubMed Abstract | Crossref Full Text | Google Scholar

24. Pegram, C, Diaz-Ordaz, K, Brodbelt, DC, Chang, YM, Hall, JL, Church, DB, et al. Later-age neutering causes lower risk of early-onset urinary incontinence than early neutering-a VetCompass target trial emulation study. PLoS One. (2024) 19:e0305526. doi: 10.1371/journal.pone.0305526

PubMed Abstract | Crossref Full Text | Google Scholar

25. O’Rourke, A, Haydock, R, Butterwick, RF, German, AJ, Carson, A, Lyle, S, et al. Estimated activity levels in dogs at population scale with linear and causal modeling. Front Vet Sci. (2025) 12:1572794. doi: 10.3389/fvets.2025.1572794

PubMed Abstract | Crossref Full Text | Google Scholar

26. Gelman, A, Carlin, JB, Stern, HS, Dunson, DB, Vehtari, A, and Rubin, DA. Bayesian data analysis. 3rd ed. London, UK: Chapman and Hall/CRC (2013).

Google Scholar

27. Spiegelhalter, DJ, Abrams, KR, and Myles, JP. In: Bayesian approaches to clinical trials and health-care evaluation. New York: John Wiley and Sons Ltd. (2004).

Google Scholar

28. R Core Team. A language and environment for statistical computing. R foundation for statistical computing. Vienna, Austria: R Core Team (2025).

Google Scholar

29. Brooks, S, Gelman, A, Jones, GL, and Meng, XL. Handbook of Markov chain Monte Carlo. Boca Raton, FL: Chapman & Hall/CRC (2011).

Google Scholar

30. Dendukuri, N, and Joseph, L. Bayesian approaches to modeling the conditional dependence between multiple diagnostic tests. Biometrics. (2001) 57:158–67. doi: 10.1111/j.0006-341x.2001.00158.x

Crossref Full Text | Google Scholar

31. Branscum, AJ, Gardner, IA, and Johnson, WO. Estimation of diagnostic-test sensitivity and specificity through Bayesian modeling. Prev Vet Med. (2005) 68:145–63. doi: 10.1016/j.prevetmed.2004.12.005

Crossref Full Text | Google Scholar

32. Berry, DA. Bayesian clinical trials. Nat Rev Drug Discov. (2006) 5:27–36. doi: 10.1038/nrd1927

PubMed Abstract | Crossref Full Text | Google Scholar

33. Yang, DA, Xiao, X, Jiang, P, Pfeiffer, DU, and Laven, RA. Keeping continuous diagnostic data continuous: application of Bayesian latent class models in veterinary research. Prev Vet Med. (2022) 201:105596. doi: 10.1016/j.prevetmed.2022.105596

PubMed Abstract | Crossref Full Text | Google Scholar

34. Woodward, AP. Bayesian estimation in veterinary pharmacology: a conceptual and practical introduction. J Vet Pharmacol Ther. (2024) 47:322–52. doi: 10.1111/jvp.13433

PubMed Abstract | Crossref Full Text | Google Scholar

35. Cummings, CO, Mitchell, MA, Perry, SM, Fleissner, N, Mayer, J, Lennox, AM, et al. Bayesian decision analysis: an underutilized tool in veterinary medicine. Animals (Basel). (2022) 12:3414. doi: 10.3390/ani12233414

PubMed Abstract | Crossref Full Text | Google Scholar

36. German, AJ, Titcomb, JM, Holden, SL, Queau, Y, Morris, PJ, and Biourge, V. Cohort study of the success of controlled weight loss programs for obese dogs. J Vet Intern Med. (2015) 29:1547–55. doi: 10.1111/jvim.13629

PubMed Abstract | Crossref Full Text | Google Scholar

37. Bryan, HM, Adams, AG, Invik, RM, Wynne-Edwards, KE, and Smits, JE. Hair as a meaningful measure of baseline cortisol levels over time in dogs. J Am Assoc Lab Anim Sci. (2013) 52:189–96. https://aalas.kglmeridian.com/view/journals/72010024/52/2/article-p189.xml

Google Scholar

38. Gao, W, Stalder, T, Foley, P, Rauh, M, Deng, H, and Kirschbaum, C. Quantitative analysis of steroid hormones in human hair using a column-switching LC-APCI-MS/MS assay. J Chromatogr B Analyt Technol Biomed Life Sci. (2013) 928:1–8. doi: 10.1016/j.jchromb.2013.03.008

PubMed Abstract | Crossref Full Text | Google Scholar

39. Raffan, E, Holden, SL, Cullingham, F, Hackett, RM, Rawlings, JM, and German, AJ. Standardized positioning is essential for precise determination of body composition using dual-energy x-ray absorptiometry in dogs. J Nutr. (2006) 136:1976S–8S. doi: 10.1093/jn/136.7.1976S

PubMed Abstract | Crossref Full Text | Google Scholar

40. German, AJ, Holden, SL, Bissot, T, Hackett, RM, and Biourge, V. Dietary energy restriction and successful weight loss in obese client-owned dogs. J Vet Intern Med. (2007) 21:1174–80. doi: 10.1892/06-280.1

PubMed Abstract | Crossref Full Text | Google Scholar

41. Broome, HAO, Woods-Lee, GRT, Flanagan, J, Biourge, V, and German, AJ. Weight loss outcomes are generally worse for dogs and cats with class II obesity, defined as > 40% overweight. Sci Rep. (2023) 13:22958. doi: 10.1038/s41598-023-50197-y

PubMed Abstract | Crossref Full Text | Google Scholar

42. van der Laan, JE, Vinke, CM, and Arndt, SS. Evaluation of hair cortisol as an indicator of long-term stress responses in dogs in an animal shelter and after subsequent adoption. Sci Rep. (2022) 12:5117. doi: 10.1038/s41598-022-09140-w

PubMed Abstract | Crossref Full Text | Google Scholar

43. Wojtaś, J, Garbiec, A, Karpiński, M, Skowronek, P, and Strachecka, A. Are hair cortisol levels of humans, cats, and dogs from the same household correlated? Animals (Basel). (2022) 12:1472. doi: 10.3390/ani12111472

PubMed Abstract | Crossref Full Text | Google Scholar

44. Stuart Marques, V, Calesso, JR, de Carvalho, OV, and da Costa-Val Bicalho, AP. Hair cortisol concentration, disease severity and quality of life in dogs with atopic dermatitis during lokivetmab therapy. Vet Dermatol. (2023) 34:339–47. doi: 10.1111/vde.13151

PubMed Abstract | Crossref Full Text | Google Scholar

45. Vernick, J, Martin, C, Montelpare, W, Dunham, AE, and Overall, KL. Understanding the influence of early-life stressors on social interaction, telomere length, and hair cortisol concentration in homeless kittens. Animals (Basel). (2025) 15:446. doi: 10.3390/ani15030446

PubMed Abstract | Crossref Full Text | Google Scholar

46. McElreath, RE. Statistical rethinking: A Bayesian course with examples in R and Stan. 2nd ed. Boca Raton, FL: CRC Press (2020).

Google Scholar

47. Pearl, J. Causal diagrams for empirical research. Biometrika. (1995) 82:669–88. doi: 10.2307/2337329

Crossref Full Text | Google Scholar

48. Greenland, S, Pearl, J, and Robins, JM. Causal diagrams for epidemiologic research. Epidemiology. (1999) 10:37–48. doi: 10.1097/00001648-199901000-00008

PubMed Abstract | Crossref Full Text | Google Scholar

49. Hernán, MA, and Robins, JM. Causal inference: what if. Boca Raton, FL: Chapman & Hall/CRC (2025).

Google Scholar

50. Textor, J. DAGitty — draw and analyze causal diagrams. (2025). Available online at: https://www.dagitty.net (Accessed July 31, 2025).

Google Scholar

51. Cinelli, C, Forney, A, and Pearl, J. A crash course in good and bad controls. Sociol Methods Res. (2022) 2022:1071–104. doi: 10.1177/00491241221099552

Crossref Full Text | Google Scholar

52. Kruschke, JK. Bayesian analysis reporting guidelines. Nat Hum Behav. (2021) 5:1282–91. doi: 10.1038/s41562-021-01177-7

PubMed Abstract | Crossref Full Text | Google Scholar

53. Heymans, MW, and Twisk, JWR. Handling missing data in clinical research. J Clin Epidemiol. (2022) 151:185–8. doi: 10.1016/j.jclinepi.2022.08.016

PubMed Abstract | Crossref Full Text | Google Scholar

54. Ryu, C. Dlookr: tools for data diagnosis, exploration, transformation. R package version 0.6.3. (2024). Available online at: https://CRAN.R-project.org/package=dlookr (Accessed July 31, 2025).

Google Scholar

55. Wickham, H, François, R, Henry, L, Müller, K, and Vaughan, D. _dplyr: a grammar of data Manipulation_. R package version 1.1.4. (2023). Available online at: https://CRAN.R-project.org/package=dplyr (Accessed July 31, 2025).

Google Scholar

56. Wickham, H. ggplot2: Elegant graphics for aata analysis. New York: Springer-Verlag (2016).

Google Scholar

57. Revelle, W (2025). Psych: procedures for psychological, psychometric, and personality research. Northwestern University, Evanston, Illinois. R package version 2.5.3. Available online at: https://CRAN.R-project.org/package=psych (Accessed July 31, 2025).

Google Scholar

58. Wickham, H. Reshaping data with the reshape package. J Stat Softw. (2007) 21:12.

Google Scholar

59. Wickham, H, and Bryan, J (2025). Readxl: read excel files. R package version 1.4.5. Available online at: https://CRAN.R-project.org/package=readxl (Accessed July 31, 2025).

Google Scholar

60. Wickham, H, Averick, M, Bryan, J, Chang, W, McGowan, LD, François, R, et al. Welcome to the tidyverse. J Open Source Softw. (2019) 4:1686. doi: 10.21105/joss.01686

Crossref Full Text | Google Scholar

61. Daniel Adler, S, Kelly, T, Elliott, T, and Adamson, J. Vioplot: violin plot. R package, version 0.1. (2025). Available online at: https://github.com/TomKellyGenetics/vioplot.

Google Scholar

62. Huntington-Klein, N. _vtable: variable table for variable Documentation_. R package version 1.4.8. (2024). Available online at: https://CRAN.R-project.org/package=vtable (Accessed July 31, 2025).

Google Scholar

63. Bürkner, PC. Advanced Bayesian multilevel modeling with the R package brms. R J. (2018) 10:395–411. doi: 10.32614/RJ-2018-017

Crossref Full Text | Google Scholar

64. Gabry, J, Simpson, D, Vehtari, A, Betancourt, M, and Gelman, A. Visualization in Bayesian workflow. J R Stat Soc Ser A Stat Soc. (2019) 182:389–402. doi: 10.1111/rssa.12378

Crossref Full Text | Google Scholar

65. Makowski, D, Ben-Shachar, M, and Lüdecke, D. BayestestR: describing effects and their uncertainty, existence and significance within the Bayesian framework. J Open Source Softw. (2019) 4:1541. doi: 10.21105/joss.01541

Crossref Full Text | Google Scholar

66. Vehtari, A, Gelman, A, and Gabry, J. Practical bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat Comput. (2017) 27:1413–32. doi: 10.1007/s11222-016-9696-4

Crossref Full Text | Google Scholar

67. Kallioinen, N, Paananen, T, Bürkner, P-C, and Vehtari, A. Detecting and diagnosing prior and likelihood sensitivity with power-scaling. Stat Comput. (2023) 34:57. doi: 10.1007/s11222-023-10366-5

Crossref Full Text | Google Scholar

68. McElreath, R. Rethinking: statistical rethinking book package. R package version 2.42, commit ac1b3b2cda83f3e14096e2d997a6e30ad109eeee. (2024). Available online at: https://github.com/rmcelreath/rethinking (Accessed July 31, 2025).

Google Scholar

69. Gabry, J, Goodrich, B, Lysy, M, and Johnson, A. Rstantools: tools for developing R packages interfacing with 'Stan'. R package version 2.4.0. 2024. (2024). Available online at: https://CRAN.R-project.org/package=rstantools (Accessed July 31, 2025).

Google Scholar

70. Ashour, SK, and Abdel-Hameed, MA. Approximate skew normal distribution. J Adv Res. (2010) 1:341–50. doi: 10.1016/j.jare.2010.06.004

Crossref Full Text | Google Scholar

71. Stan Development Team. Stan modeling language: user’s guide and reference manual. (2025). Available online at: https://mc-stan.org/docs/stan-users-guide/ (Accessed July 31, 2025).

Google Scholar

72. Abraham, SB, Rubino, D, Sinaii, N, Ramsey, S, and Nieman, LK. Cortisol, obesity, and the metabolic syndrome: a cross-sectional study of obese subjects and review of the literature. Obesity (Silver Spring). (2013) 21:E105–17. doi: 10.1002/oby.20083

PubMed Abstract | Crossref Full Text | Google Scholar

73. Van der Valk, ES, Savas, M, and van Rossum, EFC. Stress and obesity: are there more susceptible individuals? Curr Obes Rep. (2018) 7:193–203. doi: 10.1007/s13679-018-0306-y

PubMed Abstract | Crossref Full Text | Google Scholar

74. Rosmond, R. Role of stress in the pathogenesis of the metabolic syndrome. Psychoneuroendocrinology. (2005) 30:1–10. doi: 10.1016/j.psyneuen.2004.05.007

PubMed Abstract | Crossref Full Text | Google Scholar

75. Cooper, MS. Glucocorticoids in bone and joint disease: the good, the bad and the uncertain. Clin Med (Lond). (2012) 12:261–5. doi: 10.7861/clinmedicine.12-3-261

PubMed Abstract | Crossref Full Text | Google Scholar

76. Carlesso, LC, Sturgeon, JA, and Zautra, AJ. Exploring the relationship between disease-related pain and cortisol levels in women with osteoarthritis. Osteoarthr Cartil. (2016) 24:2048–54. doi: 10.1016/j.joca.2016.06.018

PubMed Abstract | Crossref Full Text | Google Scholar

77. Menger, MM, Histing, T, Laschke, MW, Ehnert, S, Viergutz, T, and Fontana, J. Cortisol stress response after musculoskeletal surgery: a narrative review. EFORT Open Rev. (2025) 10:186–92. doi: 10.1530/EOR-2024-0126

Crossref Full Text | Google Scholar

78. Bennett, A, and Hayssen, V. Measuring cortisol in hair and saliva from dogs: coat color and pigment differences. Domest Anim Endocrinol. (2010) 39:171–80. doi: 10.18637/jss.v021.i12

Crossref Full Text | Google Scholar

79. Dettenborn, L, Tietze, A, Kirschbaum, C, and Stalder, T. The assessment of cortisol in human hair: associations with sociodemographic variables and potential confounders. Stress. (2012) 15:578–88. doi: 10.3109/10253890.2012.654479

PubMed Abstract | Crossref Full Text | Google Scholar

80. Flanagan, J, Kocevar, G, Reinert, B, Morrison, J, and German, AJ. Prior obesity as a risk factor for canine diabetes mellitus: a pair-matched 1:3 case-control study. Proceedings of the 2025 ACVIM Forum, Louisville, USA (2025).

Google Scholar

81. Burnard, C, Ralph, C, Hynd, P, Edwards, JH, and Tilbrook, A. Hair cortisol and its potential value as a physiological measure of stress response in human and non-human animals. Animal production. Science. (2015) 57:401–14. doi: 10.1071/AN15622

Crossref Full Text | Google Scholar

82. Keckeis, K, Lepschy, M, Schöpper, H, Moser, L, Troxler, J, and Palme, R. Hair cortisol: a parameter of chronic stress? Insights from a radiometabolism study in guinea pigs. J Comp Physiol B. (2012) 182:985–96. doi: 10.1007/s00360-012-0674-7

Crossref Full Text | Google Scholar