Cow- and herd-level risk factors for lameness in partly housed pasture-based dairy cows

Lameness in dairy cows is a major animal welfare concern and has substantial economic impact through reduced production and fertility. Previous risk factor analyses have focused on housed systems, rather than those where cows were grazed for the majority of the year and housed only for the winter period. Therefore, the aim of this observational study was to identify a robust set of cow-level and herd-level risk factors for lameness in a pasture-based system, based on predictors from the housing and grazing periods. Ninety-nine farms were visited during the grazing period (April 2019–September 2019), and 85 farms were revisited during the housing period (October 2019–February 2020). At each visit, all lactating cows were scored for lameness (0 = good mobility, 1 = imperfect mobility, 2 = impaired mobility, 3 = severely impaired mobility), and potential herd-level risk factors were recorded through questionnaires and infrastructure measurements. Routine cow-level management data were also collected. Important risk factors for lameness were derived though triangulation of results from elastic net regression, and from logistic regression model selection using modified Bayesian information criterion. Both selection methods were implemented using bootstrapping. This novel approach has not previously been used in a cow-level or herd-level risk factor analysis in dairy cows, to the authors’ knowledge. The binary outcome variable was lameness status, whereby cows with a lameness score of 0 or 1 were classed as non-lame and cows with a score of 2 or 3 were classed as lame. Cow-level risk factors for increased lameness prevalence were age and genetic predicted transmitting ability for lameness. Herd-level risk factors included farm and herd size, stones in paddock gateways, slats on cow tracks near the collecting yard, a sharper turn at the parlor exit, presence of digital dermatitis on the farm, and the farmers’ perception of whether lameness was a problem on the farm. This large-scale study identified the most important associations between risk factors and lameness, based on the entire year (grazing and housing periods), providing a focus for future randomized clinical trials.


INTRODUCTION
Lameness is a debilitating problem in the dairy sector, representing a major welfare challenge and negatively impacting the economic sustainability of the industry (Huxley, 2012). Lameness is a painful condition that can lead to behavioral changes in dairy cows, including increased lying and decreased feeding time (Galindo and Broom, 2002). Lameness is commonly considered to be one of the 3 most costly diseases in dairy herds (Bruijnis et al., 2010), with financial losses attributable to decreased milk production and fertility as well as treatment and culling.
Risk factor studies are critical to identifying associations between potential risk factors and lameness, thus creating an important foundation for future intervention studies. Risk factors for freestall-housed cows included increased time away from the pen, decreased cow comfort, tiestall brisket boards, and no routine trimming (Espejo and Endres, 2007). In a fully housed Canadian system, herd-level risk factors included small herd sizes, slippery flooring, and reduced lying surface comfort, whereas cow-level risk factors included high parity, low BCS, and the presence of hock injuries and overgrown claws (Solano et al., 2015). Additional risk factors for cows with no pasture access include feed rail and alley design and water trough design (Sarjokari et al., 2013). Risk factors associated with lameness for cows that are grazed year-round differ from those that are fully housed. In a large-scale study in New Zealand, cow track maintenance and stockman behavior when moving animals on cow tracks were the most prominent risk factors for lameness (Chesterton et al., 1989). In an Australian pasture-based system, risk factors included rainfall levels, collecting yard stocking rate and concrete smoothness, feed-pad stocking rate, and rough handling of cattle on cow tracks (Ranjbar et al., 2016). Lameness incidence in Brazil was greater for cows that were moved faster along the cow tracks, Holstein Friesian cows compared with Jersey cows, cows with hoof abnormalities, and cows with higher parity and lower BCS (Bran et al., 2018).

Cow-and herd-level risk factors for lameness in partly
Many dairy production systems are neither fully housed nor involve year-round grazing. For example, Irish dairy farms are almost entirely spring calving and pasture based (Dillon et al., 1995); however, cows also spend approximately 4.5 mo/yr in housed facilities before calving (Dillon et al., 2019). Cows in these hybrid systems may therefore be exposed to risk factors for both systems, which may alter the relative importance of each factor. For cows in this system type, white line disease and sole hemorrhages have been reported as the most common causes of lameness (Somers and O'Grady, 2015). This is similar to cows in fully grazed systems, such as New Zealand, where noninfectious lesions were most prevalent (Chesterton et al., 2008). In contrast, cows in fully housed systems tended to have a higher prevalence of infectious lesions, such as digital dermatitis (Solano et al., 2016).
Only limited research has investigated risk factors for lameness in a part-grazed, part-housed system. Doherty et al. (2014) derived a list of potential risk factors from previous research and established how common they were in Irish herds. Somers et al. (2019) also reported cow-level risk factors as part of the same study and included higher parity, BCS loss postpartum, and lower BCS at calving. O'Connor et al. (2020) reported that herd-level factors included footbath use and holding cows in the collecting yard until milking was complete, with cow-level risk factors including stage of lactation (>120 DIM), PTA for lameness, and BCS. However, O'Connor et al. (2020) focused on the grazing period only and did not consider the housing period or evaluate directly measured farm infrastructural features as potential risk factors for lameness. Risk factors for lameness in England and Wales, where cows were out to pasture full-time in the summer, included routine trimming, use of automatic scrapers, passageway widths <3 m, and stall curb heights <15 cm (Barker et al., 2007). However, compared with the farms studied by Barker et al. (2007), Irish farms tend to be less intensive, with lower-yielding and smaller cows and a longer grazing season. Moreover, previous research involving partgrazed, part-housed systems was undertaken before or the year of the abolition of European Union milk quotas in 2015 (Barker et al., 2007;Somers et al., 2019;O'Connor et al., 2020). Therefore, opportunities have arisen for dairy farmers to undergo expansion since then (Ramsbottom et al., 2020), and potential risk factors may be altered as a result. In addition, sample sizes in these studies ranged from 10 to 49 herds; a larger-scale study would provide a more representative sample and allow a smaller effect size to be detected.
The aim of this study was to identify a robust set of the most important cow-and herd-level risk factors for lameness in a pasture-based system where cows are also housed for part of the year, using a large number of potential predictors from the grazing and housing periods. Identifying associations between risk factors and lameness will contribute to lameness prevention and deliver a focus for future intervention studies.

MATERIALS AND METHODS
Ethical approval was granted by the Teagasc Animal Ethics Committee (Cork, Ireland) before the commencement of the study (review number: TAEC202-2018). All animal measurements were carried out in compliance with the European Union (Protection of Animals Used for Scientific Purposes) Regulation 2012 (S.I. 543 no. of 2012) and the European Directive 2010/63/EU. The study involved 2 visits: one during the grazing period (April 2019-September 2019) and one during the housing period (October 2019-February 2020). The median difference between the 2 visits was 168 d [interquartile range (IQR) = 127-217], ranging from 65 to 262 d. This study was part of a larger study assessing dairy cow welfare in pasture-based systems (Crossley et al., 2021).

Farm Selection
Before recruitment of farms, selection criteria were determined to ensure that study farms represented the predominant dairy production system in Ireland; pasture based, nonorganic, and spring calving. Herds recruited had a target of ≥30 and ≤250 cows, which accounts for 95% of farms that meet the selection criteria described. Herds enrolled were registered with the Irish Cattle Breeding Federation (ICBF; Bandon, Co. Cork, Ireland); the database for all Irish-born dairy and beef cattle. Herds recruited were located within 2 h of Teagasc Moorepark for practicality reasons, and were within the main dairy farming counties of Ireland (Cork, Tipperary, Limerick, Kerry, Kilkenny, Browne et al.: LAMENESS IN PASTURE-BASED DAIRY COWS Waterford, and Wexford); 69% of all dairy cows in the country were located in these 7 counties (ICBF, 2018).
To determine the number of farms required to detect a risk factor for lameness, a simulation-based power study was performed. Multiple different scenarios were evaluated; 100 herds of 100 cows produced an estimated 93% power to detect a risk factor with a relative risk of 1.4, and 62% power to detect a risk factor with a relative risk of 1.25. A target of 100 farms was therefore deemed to be an adequate number of farms to visit.
From a list of herds provided by ICBF, 518 farms were randomly selected using SAS version 9.4 (SAS Institute Inc.), and farmers were contacted via letter or telephone to invite them to participate in the study. In total, 131 farmers responded (response rate of 25%), and 102 of these farmers were willing to participate and were deemed suitable for the study. All 102 herds were visited during the grazing period (99 farms included in statistical analysis), and 87 farms were revisited during the housing period (85 farms included in statistical analysis).

Data Collection
Details on farm management practices and facilities were collected via questionnaires and on-farm infrastructure measurements.

Farmer Questionnaire
A questionnaire was conducted with the farmer at both the first and second visit; questionnaires can be viewed as supplemental material (Browne, 2021). The questionnaire was split between the 2 visits to ensure it was not too time consuming for the farmer. Both questionnaires gathered information on the grazing and housing periods. The questionnaire at the first visit gathered information on farm background and management, cow track maintenance and grazing practices, milking practices, and lameness prevention (including routine trimming and foot bathing), detection, and treatment methods. The second questionnaire focused on housing characteristics and management, nutrition, producer demographics, and the farmers' perception of hoof health on the farm.

Infrastructure Measurements
Infrastructure measurements were taken via direct observation for the milking facilities, cow tracks, and housing facilities. Categorical scales used as part of the infrastructure measurements can be viewed as supplemental material (Browne, 2021).

Milking Facilities
Collecting yard stocking rate, presence of a slope, entrance widths, presence of a backing gate, and flooring type were recorded. The milking parlor type, size, and flooring were also recorded. At the parlor entrance and exit, the floor slipperiness (de Vries et al., 2015) and the presence of steps, slopes, sharp turns, narrow doors, and obstructions were noted. The flooring type at the parlor exit was recorded, as was the distance from the first milking unit to the end wall of the parlor, to determine the space cows had to turn after milking. The presence, type, and length of footbaths were also included in this section.

Cow Tracks
Due to time constraints, it was not possible to collect data on every cow track on each farm. Therefore, measurements were taken on the cow track in use on the day of the first visit; at the estimated halfway point between the collecting yard and the paddock, at the end point of the cow track, and at the paddock gateway. At all 3 locations the width, surface material, surface condition, presence of loose stones, and presence of a drainage ditch were recorded. The presence of loose stones was measured by placing a quadrat (0.5 m × 0.5 m), divided into 25 smaller squares, in the center of the cow track. The number of quadrat squares containing at least one loose stone was recorded. In addition, the cow track slope and camber (measured using a spirit level), the verge width, and the presence of deep wheel tracks, water erosion, and a clear channel in the road surface, suggesting a single-file path made by cows, were recorded at the end point and the halfway point.
Measurements were also taken in the segment between the collecting yard entrance and 50 m from the collecting yard along all cow tracks utilized; this was to obtain information on cow track characteristics in areas that were most regularly used by cows. At 50 m from the collecting yard, the cow track width, verge widths, and presence of loose stones were measured. The surface material, surface condition, and gradient of the steepest slope within the first 50-m segment from the collecting yard were also recorded, as well as the presence of a drainage ditch, visible slope, consistent width, sharp turns, and a single-file path made by cows.

Housing Facilities
The presence or absence of loose housing (straw yards and slatted pens) and stall housing on each farm was recorded. Housing measurements were taken in each pen that housed dairy cows. Loose housing and stall housing measurements included number of cows present at the time of the visit, accessible feed barrier length, passageway widths, flooring type, and presence of automatic scrapers and dead-ends. Lying area dimensions and the presence or absence of bedding were also noted for loose housing. For stall housing the number of stalls, overall stall condition (percentage of stalls in disrepair), and proportion of each stall type (e.g., cantilever) and direction (head-head, wall-facing, or passage-facing) were recorded. Bedding type, mat thickness, and stall hardness (McFarland and Graves, 1995) and cleanliness were also recorded for 5% of stalls (stalls randomly selected; minimum of 2 stalls). Additionally, presence of brisket board, curb height, total length, bed length, lunge space, diagonal length, neck rail height, and width were recorded for 5% of the 2 most common stall types (stalls randomly selected; minimum of 2 stalls per type).

Herd Lameness and Body Condition Scoring
All scorers undertook training with an experienced body condition scorer from Teagasc. Scorers also attended and passed a Register of Mobility Scorersapproved course in England, ensuring that lameness scoring was standardized and consistent. A total of 6 scorers were trained in body condition scoring and 4 scorers in lameness scoring. Using weighted kappa coefficients, inter-and intraobserver agreement scores were calculated for lameness scoring and body condition scoring. The mean lameness score (LS) interobserver agreement at the beginning of the first visit was 0.73 [standard deviation (SD) = 0.07], and the mean LS inter-and intraobserver agreements before the beginning of the second round of visits were 0.85 (SD = 0.06) and 0.77 (SD = 0.05), respectively. The mean BCS interobserver agreement at the beginning of the first visit was 0.74 (SD = 0.06), and the mean BCS inter-and intraobserver agreements before the second visit were 0.81 (SD = 0.06) and 0.87 (SD = 0.05), respectively.
Herd scoring was carried out after milking at each visit; cows were retained in a crush (chute) to enable tag number identification and body condition scoring. At both visits, the number of cows in the milking herd to assess for BCS was calculated based on herd size using the Welfare Quality sample size protocol (Welfare Quality Consortium, 2009). The cows were scored using a scale from 1 to 5, in 0.25 increments (Wildman et al., 1982), by one observer. All cows in the milking herd were subsequently individually scored for lameness as they left the crush, by a single observer using the Agricultural and Horticultural Development Board Dairy 4-point scale (Archer et al., 2010).

Herd Management Data
Cow-level data were provided by the ICBF for all herds enrolled in the study. Date of birth and date of first calving were classified into age at visit (yr) and age at first calving (mo), respectively. Days in milk on the day of the visit, calving interval (between 2018 and 2019 calving), and days until next calving were calculated based on calving dates provided. Based on the 2019 lactation, the parity, calving difficulty, whether the cow had twins or a single calf, average SCC, 305-d milk recording prediction, and dry-off date were provided for each cow. Breeds were classified into Holstein Friesians, other purebreds (excluding Holstein Friesians), and crossbreeds. Purebreds were defined as cows that were ≥87.5% of a single breed. The 2019 Economic Breeding Index, maintenance subindex, and health subindex values were extracted for each cow; explanations of these indices can be found in Berry et al. (2007). The lameness trait within the health subindex, in the form of a PTA, was also provided. In terms of lameness, a positive PTA indicates that the progeny are more likely to become lame than the base population (Berry et al., 2007).

Statistical Analysis
All data cleaning, pre-processing of data, descriptive statistics, and statistical modeling were executed in R software version 3.3.1 (R Core Team).

Data Cleaning
A total of 22,164 LS were recorded across 102 farms. Three farms, comprising 262 LS observations, were excluded from the data set due to robotic milking (1 farm) or once-a-day milking (2 farms). A further 1,694 LS observations were removed due to wrongly recorded tag numbers, accidental scoring of pre-calving heifers, and scoring of non-spring-calving cows.

Pre-Processing
Before statistical analysis, all housing predictors were weighted by the number of cows in each pen, to account for varying number of cows being subjected to the conditions of each pen. Continuous cow-level variables with missing values were split into quartiles, and an additional category was made for both continuous and categorical variables, to represent missing data points (<1% of data set). Nonparametric methods based on random forest algorithms were employed to impute missing values (3.2% of data set) from the surveys and on-farm measurements, using the missForest package (Stekhoven, 2013). Thirty-five predictors with near-zero variance were removed (Kuhn, 2020), leaving a final data set consisting of 197 predictors (cow-level predictors = 16; herd-level predictors = 181). These predictors can be viewed as supplemental material (Browne, 2021). Continuous predictors were subsequently centered and scaled (to SD units relative to overall mean). Each cow was assigned a lameness outcome at each visit: LS of 0 or 1 was classified as non-lame, and a score of 2 or 3 was classified as lame.

Variable Selection Models
Triangulation (Lawlor et al., 2017;Lima et al., 2021) of results from elastic net regression (Enet), a form of regularized logistic regression, and logistic regression using modified Bayesian information criterion (mBIC) was used to establish important risk factors for lameness. These methods were chosen due to the large number of predictors and the need to avoid overfitting. The outcome variable was lameness status (0 = not lame, 1 = lame); lameness scores from both the grazing and housing visits were included in the models. All covariates described previously were offered to each model.

Elastic Net Regression
Elastic net regression combines the ridge penalty (penalizing the sum of squared coefficients) with the lasso penalty (penalizing the sum of coefficients). The elastic net penalty term is shown in Equation [1]: where λ is a model tuning parameter providing coefficient penalization; α is the mixing parameter to determine the proportion penalty applied as ridge or lasso, where α = 0 represents a full ridge model and α = 1 represents a full lasso model; j represents a predictor variable and P represents the total number of predictors; β represents the sum of coefficients. Elastic net regression was performed using the packages caret (Kuhn, 2020) and glmnet (Friedman et al., 2010). An Enet model was fitted using a large tuning grid of α values (α = 0, 0.2, 0.4, 0.6, 0.8, 1) and λ values (λ = 0.0001, 0.001, 0.003, 0.004, 0.01, 0.015, 0.02, 0.03, 0.04, 0.05, 0.1). Five-fold cross validation with 10 repeats was used to evaluate model performance and select the best-performing model based on accuracy.

Selection Using Modified Bayesian Information Criterion
Any predictors not correlated with the outcome variable (Pearson correlation test, P > 0.3) were removed, and stepwise logistic regression model selection based on minimizing mBIC was performed using the bigstep package (Szulc, 2019). The model was fitted to best balance the penalty term against a measure of model fit. The mBIC penalty term can be described as follows (Equation [2]): where k i represents the number of predictors in the ith model, n is the number of observations, and p represents the probability that a predictor, chosen at random, influences the outcome variable (lameness status).

Bootstrapping
Bootstrapping was used for both the Enet and mBIC model selection processes. One thousand bootstrap repeats were performed for each model type; this was deemed sufficient to obtain an accurate 95% bootstrap percentile confidence interval (Efron and Tibshirani, 1993). For each bootstrap repeat, the coefficient for each predictor was returned, and the mean of the nonzero coefficients and the 95% bootstrap confidence interval for each predictor was calculated. Coefficient means were subsequently unstandardized by dividing by the SD, and odds ratios (OR) calculated using these values.

Stability Selection and Model Triangulation
A stability value was calculated for each predictor for each model selection method (elastic net regression and selection based on mBIC), defined as the proportion of bootstrap repeats in which the coefficient for that predictor was nonzero. A nonzero coefficient implied that the predictor was selected in the model. A bootstrap P-value was also determined for each predictor based on the distribution of nonzero coefficients. The P-value was calculated as the proportion of coefficients on the minority side of zero.
Drawing on the principles of stability selection, for which it is known that variables with the highest stability values are least likely to be false positives (Lima et al., 2021;Meinshausen and Bühlmann, 2010), and triangulation, for which it is accepted that use of multiple analyses reduces bias in results (Lawlor et al., 2017), final model selection was based on high-stability variables that occurred in both model types. Predictor variables that had a bootstrap P-value of <0.05 and were ranked in the top 24 by stability (number of predictors that had a stability of >80% in the Enet model, a previously established technique: Lima et al., 2020) for each method were deemed likely to have important associations with lameness. The final subset of results was not found to be sensitive with the arbitrary choice of selecting predictors ranked in the top 24 by stability. An identical subset of predictors was found if selection was based on the top 30 predictors ranked by stability.

Potential Clustering Effect
The effect of accounting for clustering at herd level and cow level were evaluated by estimating parameters for random effects logistic regression models using Markov chain Monte Carlo via the brms package (Bürkner, 2017). One model included a random effect representing herd, and a second model included random effects terms at cow and herd level. A subset of covariates was included in the logistic regression models based on those selected in both the Enet and the mBIC models. Coefficients from each model were assessed to ensure that direction of association was the same as (and effect size similar to) the results from triangulation of the Enet and mBIC models.

Cow Characteristics and Lameness Prevalence
The median age across all cows scored was 5 yr (IQR 3-7) with a median parity of 3 (IQR 2-5). The median 305-d milk yield was 6,638 kg per cow (IQR 5,750-7,597) with a median calving interval of 369 d . The median BCS during the grazing visit and the housing visit were 3 (IQR 3-3.25) and 3.25 (IQR 3-3.5), respectively. Of all cows scored, 51% were Holstein Friesian, 28% were crossbreeds, and 21% were other purebreds. The final data set consisted of 20,208 LS recorded across 99 farms. Cow-level lameness prevalence (LS2 and LS3) was 9.3% during the grazing period and 8.9% during the housing period. Lameness prevalence across farms ranged from 0.9% to 31.4% during the grazing period and from 0% to 28.0% during the housing period. Figure 1 shows the stability and bootstrap P-value for each predictor in both Enet and mBIC models, illustrating variables selected in the triangulation pro-cess. Twenty-four predictors were selected in the final Enet and final mBIC models. Of these predictors, 11 were selected in both models and therefore represented a robust set of risk factors for lameness. Figure 2 shows the standardized mean coefficient and 95% confidence intervals for each predictor that was selected in both models (for comparison of effect size). Table 1 reports full results for predictors selected in both the Enet and mBIC models. Random effects logistic regression models suggested that accounting for the clustering effect of (1) herd and (2) cow nested within herd did not substantially influence the results from the Enet and mBIC models.

Cow-Level Risk Factors
Age had the largest standardized effect size of all cow-and herd-level predictors (based on the average of the standardized mean coefficients from the Enet and mBIC models); as age increased by 1 yr, the odds of a cow being lame increased by approximately 20% (Enet OR = 1.19; mBIC OR = 1.21; mean OR = 1.20). A positive lameness PTA increased the odds of lameness by approximately 37.5% (Enet OR = 1.14; mBIC OR = 1.61) compared with those with a negative PTA.

Herd-Level Risk Factors
Five herd-level factors were associated with an increased risk of lameness. In both the Enet and mBIC models, "farmers who considered lameness to be a problem in their herd" had the largest standardized effect size of all herd-level predictors. When farmers considered lameness to be a problem in their herd, odds of lameness increased by approximately 47% (Enet OR = 1.17; mBIC OR = 1.77) compared with when farmers did not consider lameness to be a problem. When ≥10% of the herd had been treated for lameness in the year before the study, the odds of lameness were increased by approximately 27% (Enet OR = 1.08; mBIC OR = 1.46) compared with those herds where <10% were treated. Additionally, when >5% of the herd had digital dermatitis during the current lactation, according to the farmer, the odds of lameness were increased by approximately 30% (Enet OR = 1.08; mBIC OR = 1.52) compared with a herd with ≤5%. A 10% increase in the proportion of slats in the first 50 m of cow tracks following the collecting yard increased the odds of lameness by approximately 6.5% (Enet OR = 1.04; mBIC OR = 1.09). Also, a 10% increase in the percentage of the gateway surface material that was stones increased the odds of lameness by approximately 7% (Enet OR = 1.03; mBIC OR = 1.11).
Four herd-level predictors reduced the risk of lameness. As herd size increased by 100 cows, the odds of lameness decreased by approximately 23% (Enet OR = 0.90; mBIC OR = 0.64). Similarly, as the grazing platform size increased by 100 ha, the odds of lameness decreased by approximately 45% (Enet OR = 0.74; mBIC OR = 0.36). Herds with no digital dermatitis cases during the current lactation, according to the farmer, had decreased odds of lameness of approximately 20.5% (Enet OR = 0.91; mBIC OR = 0.68) compared with a herd with >0% and ≤5%. Also, as the distance to turn after milking increased by 1 m, the odds of lameness decreased by approximately 8.5% (Enet OR = 0.97; mBIC OR = 0.86).

Predictors Selected in Individual Models but Excluded in Triangulation
A larger set of predictor variables were selected in one or the other of the individual models (Enet or mBIC), but not in both, and were therefore not reported in the previous triangulated results. Thirteen of the predictors that were selected in the final Enet model (within the top 24 ranked by stability) were not also selected in the final mBIC model. These predictors, which showed increased risk for lameness, included a higher proportion of cubicles of recommended width, higher proportion of cow track surface material measurements recorded as stones, and longer walking distance to the furthest paddock. Predictors that showed decreased risk for lameness included higher economic breeding index, genetic health and maintenance subindexes, greater days in milk, first-parity cows (i.e., no calving interval, compared with cows with a calving interval of 353 to 369 d), higher proportion of cubicles with thick mats, higher proportion of collecting yard that was grooved concrete, higher proportion of cow track measurements with a gradient >10%, higher proportion of cows tracks in the first 50 m from the collecting yard with a ditch, and higher proportion of cow tracks where the transition from concrete to other surface material was within 50 m of the collecting yard entrance.
Similarly, 13 of the predictors for lameness that were selected in the final mBIC model were not also selected in the final Enet model. Predictors that showed increased risk for lameness included high number of  Table 1. Results from the elastic net regression model (Enet) and from the logistic regression model using modified Bayesian information criterion (mBIC), ordered by mean stability rank across both models, to identify risk factors for lameness from 99 spring-calving, pasture-based herds during the grazing period (April 2019-September 2019) and in 85 of these herds during the housing period (October 2019-February 2020); the mean standardized coefficient and 95% CI, odds ratio, stability rank, and bootstrap P-value are shown for each predictor selected in both the Enet and mBIC models; the binary outcome is lameness status (0 = not lame, 1 = lame) Predictor  Odd ratios calculated from unstandardized means. For continuous predictors, a unit change of 0.1, 10, and 100, as well as the usual 1-unit change, was used to enable result interpretation.

2
A positive lameness PTA indicates that the progeny are more likely to become lame than the base population (Berry et al., 2007). 3 During the current lactation, according to the farmer. 4 In the year before the study started (2018), according to the farmer. stones in paddock gateway, dry cow cubicles cleaned once per day (compared with less than once per day), cows housed based on parity, cow tracks repaired less than once per 2 years (compared with once per year or more frequently), PTA of 0 for lameness (compared with negative PTA), mobility scoring visit carried out in May (compared with April), BCS <3 (compared with BCS 3), all cow track points measured as wide enough based on herd size (compared with farms with a combination of cow tracks measured that were wide and narrow), and herds that were routinely trimmed. Predictors that showed decreased risk for lameness included second-parity cows (compared with first-parity cows), third-parity cows (compared with first-parity cows), copper sulfate and formalin used in foot bathing routine (compared with copper sulfate only), and BCS >3.25 (compared with BCS 3).

DISCUSSION
From a cohort of approximately 200 potential cowlevel and herd-level predictors in the final model, 11 risk factors were deemed highly likely to have important associations with lameness in partly housed pasture-based dairy cows. To the authors' knowledge, this is the first time important lameness predictors have been found based on the entire year, in this particular system, and the first time this novel statistical approach (mBIC and Enet triangulation with bootstrapping) has been used to identify risk factors in dairy cows.

Cow-Level Risk Factors
In agreement with previous studies, the risk of lameness increased with age (Rowlands et al., 1985;Haskell et al., 2006). This may be explained by changes in the functional anatomy of the hoof with age, such as the degeneration of the digital cushion (Räber et al., 2004). Irreversible bone development on the distal phalanx has also been reported to increase with age, history of lameness, and previous cases of sole ulcers, sole hemorrhages, and white line disease (Newsome et al., 2017). Additionally, older cows are more likely to have a history of lameness, and previous lameness has been shown to be a major predictor of future lameness (Randall et al., 2018). In contrast to the current study, the study by Randall et al. (2018) had a longitudinal study design and therefore provides much stronger evidence for causality. However, the study included only 2 UK dairy farms and may not be comparable to Irish dairy farms, where all cows have prolonged pasture access and cows are generally lower yielding. Aging is inevitable; however, the effect of aging on lameness can be minimized through prevention of first-time lameness events, early detection of lameness, and effective treatment of lesions (Randall et al., 2018).
Cows with a positive lameness PTA exhibited a higher risk for lameness than cows with a negative PTA. Lameness PTA is a specific genetic index, in which a higher lameness PTA indicates the progeny will have a higher susceptibility to lameness (Berry et al., 2007). O 'Connor et al. (2020) reported similar findings: a positive lameness PTA compared with a negative PTA increased the odds of lameness by 41%. Similarly, the current study showed an increased odds ratio of 44%. These results add support for the lameness PTA and emphasize that genetic selection is influential for reducing lameness at cow level. The choice of bulls used for breeding may be more important as a long-term lameness reduction strategy than previously realized.

Herd-Level Risk Factors
The results of this study provide no evidence that farm expansion increases the risk of lameness in a partgrazed, part-housed system. As reported previously, a larger herd reduced lameness risk (Dippel et al., 2009;Chapinal et al., 2013). Solano et al. (2015) reported that a herd size of more than 100 cows reduced the odds of lameness by one-third, compared with a herd size of less than 100. Despite cows walking longer distances on larger pasture-based farms (Beggs et al., 2019), improved management and facilities could explain the reduced lameness prevalence. In contrast, Alban (1995) reported that lameness was positively correlated with herd size, which may be explained through more cows per staff member (Sundrum, 2015) and poorer recognition of individual cows (Dippel et al., 2009) in larger herds. Other studies have also reported that herd size was not significant in relation to lameness (Espejo and Endres, 2007;Barker et al., 2010;Beggs et al., 2019). The varied results highlight the lack of clarity regarding the association between herd size and lameness, and the interplay with other factors that influence this relationship.
Lameness risk was reduced when cows had a longer distance to turn at the parlor exit. All parlors in this study were herringbone or parallel, meaning that cows exited the parlor in single file, usually making a 90-or 180-degree turn onto a passageway to return to their pasture or pen. Similarly, Barker et al. (2010) reported that sharp turns at the parlor entrance or exit increased the risk of lameness. Similarly to the current study, the cross-sectional study design used by Barker et al. (2010) does not prove a causative relationship between sharp turns and lameness; however, it can be used to establish causal hypotheses. One commonly posited theory is that shearing forces on the hoof when cows Browne et al.: LAMENESS IN PASTURE-BASED DAIRY COWS turn sharply can lead to white line disease, potentially explaining the negative correlation between turning distance and lameness prevalence. Sharp turns may also reduce cow flow, instigating crowding and pushing of cows at the parlor exit. Rubber matting has been proven to increase friction and compressibility, in turn reducing slipping and improving mobility (Rushen and de Passillé, 2006). Therefore, introducing rubber matting where sharp turns are present at the parlor exit may be beneficial in improving cow flow and reducing claw trauma. Randomized clinical trials proving the effectiveness of this intervention are currently lacking; further research in this area is required.
Slats in the first 50 m of cow track following the collecting yard increased the risk of lameness in this study. Slatted flooring has previously been linked to increased lameness prevalence (Dippel et al., 2009) and claw health problems (Burgstaller et al., 2016) in housed cattle; however, limited information exists on the implications of slats on cow tracks. Concrete slats are more slippery compared with solid concrete flooring (Rouha-Mülleder et al., 2009), leading to a reduced pace and shortened strides (Telezhenko and Bergsten, 2005). Slatted flooring also creates uneven weight distributions across the claws, predisposing to white line disease (Hinterhofer et al., 2006). Installing rubber matting over the slats could reduce slipperiness, hoof lesions, and overall lameness prevalence (Hultgren and Bergsten, 2001;Telezhenko and Bergsten, 2005).
Stones in the gateways to pasture also presented a risk for lameness. It is hypothesized that stones penetrate the hoof horn, causing separation of the white line, and subsequently lead to an infection of the dermal tissue in more severe cases. An uneven stony surface may also result in shearing forces on the hooves. Although this study provides no evidence for causality, the association identified between stones in gateways and lameness supports these theories. Gudaj et al. (2012) reported that cows required more blocks during trimming when stones were present on cow tracks. However, in contrast to the current study, all cows in the study by Gudaj et al. (2012) were Holstein Friesian, and only cows on 14 farms, out of 25 farms visited, had access to pasture. Gateways may be more high-risk areas due to cows pushing through a narrow entrance and being unable to avoid stones. Where finances are limited, it may be beneficial to prioritize maintenance of commonly used gateways, to ensure minimal stones are present, before general cow track maintenance.
Three of the risk factors identified are subjective impressions of the farmer: the presence of digital dermatitis in the herd, the percentage of the herd treated for lameness, and farmers who consider their herd to have a lameness problem. Although these results are not entirely unexpected, they indicate that farmers in this study acknowledge lameness as an issue and can therefore work toward eliminating the disease. This is in contrast to previous studies, which have generally shown that a low proportion of farmers perceive lameness to be a problem in their herd Sadiq et al., 2019).
Based on predictors identified by both the final mBIC and the Enet model, no characteristics specifically linked to housing infrastructure and management were found to be important risk factors for lameness in a typical Irish dairy system. This emphasizes that the housing period did not seem to have a large influence on lameness, in contrast to the grazing period. Cows are only housed for approximately one-third of the year in Ireland; therefore, cows are exposed to the effects of grazing for a more prolonged period of time, and thus the grazing period appears to have the greatest influence on lameness development. However, although the grazing period was shown to have the greatest influence on lameness development, some housing features and management were selected in one or another of the Enet and mBIC models (although excluded by triangulation as not selected in both), such as cubicle mat thickness and frequency of cubicle cleaning for dry cows. Due to these variables not being selected in both the Enet and mBIC models via triangulation, it is less likely that these are generalizable for the target population; these predictors may have smaller effect sizes and may be very important on some farms and not in others.

Modeling Methods
The multifactorial nature of lameness and the need to construct a statistical model based on a large number of predictor variables would likely lead to problems with overfitting in simple regression models (Vatcheva et al., 2016); this is increasingly recognized as a potential feature in a large proportion of previous work across a range of disciplines. This is especially problematic where the sample size is small relative to the number of potential predictors; in this case, the sample of lameness scores was relatively large, but the vast majority of predictors varied only at farm level. To overcome this issue, regularized regression (Zou and Hastie, 2005) and selection using mBIC (Bogdan et al., 2008) have both been proposed for variable selection. As the ability to capture large amounts of data on-farm improves and data sets become wider, these methods will become increasingly important in statistical analysis. Using conventional methods such as stepwise selection based on Akaike's information criterion, a larger set of risk factors would likely have been identified that were false positives and likely to have inflated coefficients (Hastie et al., 2015;Lima et al., 2021). In this study, a relatively conservative analytical approach was chosen, to minimize the chances of reporting false-positive risk factors. The additional predictors included in one or the other of the 2 models represent a set of factors that can more speculatively be associated with the outcome, and it is worth noting that these would have been reported as significant predictors had a single modeling approach been chosen. The aim of this study was to identify a set of risk factors that are the most important in a pasturebased system and are most likely to be generalizable across a high proportion of similar farms.
Between-model variation was also accounted for through triangulation (Lawlor et al., 2017) of the Enet and mBIC models. Triangulation combines results from multiple statistical methods to obtain reliable results, because the bias from each model type is discounted (Lawlor et al., 2017;Lima et al., 2021). Elastic net regression has a tendency toward a higher false-positive rate and deflated coefficient values, whereas mBIC has a higher false-negative rate and inflated coefficient values (Lima et al., 2021), displaying opposing biases. The difference in effect size between the Enet and mBIC models observed for some predictors is therefore not unexpected, and it is likely that the true effect size lies in between the 2 estimates. These methods have allowed identification of a robust list of risk factors and direction of effect, and have given an indication of likely effect size.
Within-model variation was also accounted for through bootstrapping, a resampling technique for statistical inference (Dixon, 2002). Bootstrapping is beneficial to assess variable stability and coefficient distribution (Sauerbrei and Schumacher, 1992;Meinshausen and Bühlmann, 2010). To the best of the authors' knowledge, regularized regression and model selection using mBIC, with the use of bootstrapped selection stability, have not previously been used in a cow-level or herd-level risk factor analysis among dairy cows.

Study Limitations
This study may be susceptible to some bias due to farmers having the opportunity to choose whether to participate in the study. However, a selection criterion was established before recruiting participants, to ensure farms were representative of a typical Irish dairy farm. Additionally, several of the observations and measurements were slightly subjective, therefore leading to potential bias. This study has a cross-sectional design, and, as such, the associations found do not imply causation. This study design is valuable for assessing a large number of potential risk factors at once, without the logistical challenges of running multiple expensive ran-domized clinical trials. Lameness typically occurs after exposure to a risk factor; therefore, exposure to a risk factor during the end of the grazing period may result in lameness during the subsequent housing period, and, similarly, exposure to housing risk factors may result in lameness during the subsequent grazing period. This issue was acknowledged by including lameness scores and potential predictors from both the housing and grazing periods in the same model. This also allowed the most important risk factors, based on the entire year, to be identified. Findings from this study provide a base of knowledge and deliver a focus for future lameness intervention studies in Irish pasture-based systems.

CONCLUSIONS
Both cow-level and herd-level risk factors were associated with lameness in a part-grazed, part-housed system. Triangulation of bootstrapped regularized regression and logistic regression model selection based on modified Bayesian information criterion proved a robust way to identify a subset of important risk factors from a very large number of potential predictors. Cow-level risk factors included increased age and a positive PTA for lameness. Herd-level risk factors included smaller herd size and grazing platform, increased presence of digital dermatitis, presence of stones in gateways and slats on cow tracks, a tighter turn following milking, farmers who treated a higher proportion of their herd for lameness, and farmers who considered lameness to be a problem in their herd. Based on this study, farmers may benefit from a breeding program that places greater emphasis on lameness traits, taking measures to mitigate the effect of tight turns at the parlor exit and slats on the cow tracks, and removing stones from paddock gateways. Applying a package of measures across multiple herds in a randomized clinical-type trial, such as putting matting at the milking parlor exit and replacing slats on the cow tracks, might be useful for determining effective methods for decreasing lameness in Irish dairy cows.
Medicine and Science, University of Nottingham, Nottingham, UK) for statistical advice. We also acknowledge the Irish Cattle Breeding Federation (Bandon, Ireland) for their contribution. This study was funded by the Teagasc Walsh Scholarship Program (Carlow, Ireland) and Dairy Research Ireland. The authors have not stated any conflicts of interest.