The analysis was conducted at the herd level using convenience samples from \(n = 10\) cattle herds. For each herd \(i\), the infection outcome was summarized as a binomial count:

\[ Y_i \sim \text{Binomial}(n_i, p_i), \qquad Y_i = k_i \text{ infected out of } n_i \text{ tested}, \]

where \(p_i\) represents the herd-level probability of gastrointestinal nematode (GIN) infection. This formulation accounts for both the number of animals tested per herd and the proportion infected, and allows inference directly at the herd level—an appropriate unit given that management practices (e.g., stocking density and deworming) are implemented at that scale.

Primary predictors were selected based on biological relevance and prior literature linking herd-level infection pressure to environmental and management factors. The following variables were considered:

Two biologically justified interaction terms were specified to allow for environmental modification of the herd-size effect: \[ z_{\text{TA}} \times \texttt{elev}, \quad z_{\text{TA}} \times \texttt{den}. \] These interactions test whether the relationship between herd size and infection risk differs between elevation or density strata, which could reveal context-dependent effects of herd management.

Reference categories were defined so that the intercept corresponds to the expected log-odds of infection for an average-sized herd (\(z_{\text{TA}} = 0\)) located at high elevation, under high stocking density, and without deworming.

The final model () was a binomial logistic regression of the form: \[ \text{logit}(p_i) = \beta_0 + \beta_1 z_{\text{TA},i} + \beta_2 1(\text{elev} = \text{low})_i + \beta_3 1(\text{den} = \text{low})_i + \beta_4 1(\text{deworm} = \text{Yes})_i + \beta_5 [z_{\text{TA},i} \times 1(\text{elev} = \text{low})_i] + \beta_6 [z_{\text{TA},i} \times 1(\text{den} = \text{low})_i]. \]

Models were estimated by maximum likelihood using the canonical logit link. For each herd \(i\), the binomial likelihood is given by: \[ L_i(\boldsymbol{\beta}) = \binom{n_i}{k_i} p_i^{k_i} (1 - p_i)^{n_i - k_i}, \qquad p_i = \frac{1}{1 + \exp(-\eta_i)}, \quad \eta_i = \mathbf{x}_i^{\top}\boldsymbol{\beta}. \]

This structure allows infection probability to vary simultaneously with herd size and environmental conditions while preserving interpretability of the coefficients as log-odds effects.

Primary inference relied on Wald \(z\)-tests for individual coefficients using two-sided \(\alpha = 0.05\). Effect sizes were expressed as odds ratios (OR) with 95% confidence intervals: \[ \text{OR} = \exp(\hat{\beta}), \quad 95\%\ \text{CI} = \exp(\hat{\beta} \pm 1.96 \times \text{SE}_{\hat{\beta}}). \] Given that interactions were pre-specified based on biological reasoning, interpretation emphasized the simple slopes of herd size (\(z_{\text{TA}}\)) within each environmental stratum and the direction of effect modification rather than overall main effects.

Model explanatory value was assessed through an analysis of deviance comparing the full model against the intercept-only null model. The deviance difference, \(\Delta D\), was evaluated against the chi-square reference distribution using the corresponding difference in degrees of freedom (\(\Delta df\)). This test provides a global assessment of whether inclusion of the predictors and interactions significantly improves model fit.

To evaluate whether the binomial variance assumption was appropriate, the dispersion parameter was estimated as \[ \hat{\phi} = \frac{\text{Residual deviance}}{\text{Residual degrees of freedom}}. \] Values near unity indicate adequate dispersion; \(\hat{\phi} > 1\) suggests overdispersion (extra-binomial variation), whereas \(\hat{\phi} < 1\) implies underdispersion. For model , the estimated dispersion (\(\hat{\phi} = 0.69\)) indicated no evidence of overdispersion and therefore no need for quasi-binomial adjustments.

Standard influence diagnostics were examined to ensure stability of parameter estimates and identify potential high-leverage herds.
Cook’s distance (\(D_i\)), DFFITS, and DFBETAs were computed using the built-in function in R. Heuristic cutoffs used for flagging potential influence were \(D_i > 4/n\), \(|\text{DFFITS}_i| > 2\sqrt{p/n}\), and \(|\text{DFBETA}_{ij}| > 2/\sqrt{n}\), where \(p\) is the number of model parameters and \(n = 10\) herds.

To assess robustness, leave-one-herd-out refits were conducted. Each iteration removed one herd, refit the model, and compared the resulting coefficients with the full-data estimates. Stability of coefficient signs and relative magnitudes was interpreted as evidence of robustness. For predictive evaluation, the binomial log-likelihood of the held-out herd was computed at each step, and the mean held-out log-likelihood (\(-3.00\) for ) was compared to the alternative elevation specification (\(-7.18\) for ). The higher (less negative) mean log-likelihood confirmed better predictive performance for the categorical-elevation model.

Cases showing extreme predicted probabilities (approaching 0 or 1) were interpreted as evidence of quasi-separation in sparse strata. To verify that these patterns did not distort inference, a bias-reduced (Firth-type) logistic regression was fit as a sensitivity check. The direction of all coefficients was consistent with the maximum-likelihood estimates, confirming robustness of the substantive conclusions.

Predicted probabilities and 95% confidence intervals were obtained by transforming the fitted linear predictor via the inverse-logit function: \[ \hat{p}_i = \frac{1}{1 + \exp(-\hat{\eta}_i)}, \qquad \hat{\eta}_i \pm 1.96\, \text{SE}(\hat{\eta}_i) \xrightarrow{\text{inv-logit}} \text{CI on } p_i. \] Predictions were generated using from the package to summarize marginal effects and interactions.
Plots were produced for the full combination of elevation () and density () across a grid of standardized herd sizes, stratified by deworming status.
Two-way slices (herd size by elevation, herd size by density) were also presented with other factors held at reference conditions.
This approach isolates the shape and direction of interaction effects while avoiding extrapolation beyond the observed \(z_{\text{TA}}\) range.

All herds had complete data for outcome and predictors; no imputation was required.
Continuous predictors were standardized before fitting to enhance comparability and mitigate collinearity in interaction terms.
Categorical reference levels (high elevation, high density, no deworming) were chosen to represent the most prevalent and ecologically meaningful baseline.
Given the limited sample size (\(n = 10\)), the model was restricted to pre-specified main effects and two biologically justified interactions to avoid over-parameterization and inflated variance.

Because of small sample size and sparsity in certain combinations (e.g., few low-density herds and few untreated herds), interpretation emphasized the direction of effects and overlap of confidence intervals rather than precise magnitude estimates in underrepresented strata.
Consistency across multiple diagnostics (Cook’s \(D\), DFFITS, DFBETAs) and leave-one-herd-out analyses was prioritized as evidence of model stability.
Comparative predictive performance between the continuous and categorical elevation formulations further supported selection of the categorical specification () as the final model, balancing interpretability, ecological coherence, and diagnostic behavior.

All analyses were performed in (version to be specified in the final document).
Binomial logistic models were fit using the base function with .
Predicted probabilities and confidence bands were computed from the fitted linear predictor using the inverse-logit transformation.
Visualizations were generated within a grammar-of-graphics framework using , and all results are fully reproducible from the code and data described above.