Resolving design of experiments for factorial layouts with applications to fraser valley dairy farm productivity

Introduction

In 1926, R. A. Fisher published the first systematic account of factorial experimentation, arguing that testing multiple factors simultaneously is more efficient than the one-factor-at-atime approach that had dominated agricultural trials for a century [5] . His insight reshaped experimental science, and factorial designs are now standard in fields from pharmaceuticals to semiconductor manufacturing. But one old question persists: how few runs can you get away with when the number of factor levels grows beyond two or three? The resolution of a fractional factorial design describes which effects are confounded with one another [6] . A resolution III plan confounds main effects with two-factor interactions; resolution IV separates main effects but confounds pairs of two-factor interactions; and resolution V or higher keeps all main effects and two-factor interactions estimable. For 2^k factorials, the theory is mature: defining relations, word-length patterns, and minimumaberration criteria provide a complete toolkit for selecting fractions [7] . For factorials with more than two levels-say, 5^4, as in the dairy application considered here-the situation is less tidy. Latin square and Graeco-Latin square constraints can be overlaid on the design to control for nuisance variables, but finding the fraction that maximises information while obeying these constraints is generally an NP-hard combinatorial problem [8] .

Convex optimisation offers a way around the combinatorial barrier. By relaxing the integer requirement on the design measure, one can solve a continuous log-determinant problem whose solution provides an upper bound on the information achievable by any exact design of a given size [1] . Rounding heuristics then convert the continuous measure into an implementable list of runs, and the gap between the rounded design's D-efficiency and the continuous bound quantifies how much was lost in the conversion. This approach has been applied to 2-level and 3-level factorials but rarely to designs with 5 or more levels [9] . The Fraser Valley, located east of Vancouver, supports one of Canada's densest dairy regions, with about 430 farms producing over 800 million litres of milk annually [10] . Feed composition-specifically the proportions of corn silage, alfalfa hay, canola meal, and a mineral premix-is a controllable factor that directly affects milk yield, fat percentage, and somatic cell count. Testing four feed factors at five levels each would require 625 runs in a full factorial, an impractical number when each "run" means feeding a herd a particular diet for an entire milking season. This research used the convex relaxation approach to design a fractional factorial with only 50 runs (8% of the full factorial) and then evaluated its performance on data from 16 Fraser Valley farms over two seasons (2020-2022). Trimmed means handled outlying fat-percentage observations, and a hazard function tracked the time to production drop-off in individual cows. A cryptographic hash of the randomisation sequence served as a tamper-detection measure, a practical safeguard in multi-farm trials where blinding is impossible [4] .

Model Description and Parameterisation

The information matrix for a k-factor design with n_i levels per factor is M(ξ) = ∑ ξ_j f(x_j) f(x_j)ᵀ, where ξ_j is the weight assigned to design point x_j and f(x_j) is the vector of model terms evaluated at that point [1] . For a main-effectsplus-interactions model with four factors at five levels, f has 1 + 4(4) + 6(16) = 113 entries. The D-optimality criterion maximises log det M(ξ) subject to ∑ ξ_j = 1, ξ_j ≥ 0. Latin square constraints were encoded as ∑_{x: x_k = l} ξ_j = 1/n_k for each factor k and level l. The interior-point solver SDPT3, called through the CVX modelling framework in MATLAB R2023a, solved the relaxed problem in under 40 seconds on a desktop workstation. Rounding to an exact 50run design used the Fedorov exchange algorithm seeded with the 50 support points of highest ξ-weight.

Sensitivity Analysis

To assess how robust the D-efficiency was to perturbations in the number of runs, the convex relaxation was solved for design sizes ranging from 30 to 100 in steps of 10. D-efficiency rose steeply from 0.72 at 30 runs to 0.89 at 50 runs, then flattened, reaching only 0.93 at 100 runs. The marginal gain from each additional 10 runs dropped below 0.01 beyond 70 runs, confirming that 50 runs sat near the knee of the efficiency curve. Factor-level balance, measured by the maximum deviation from equal replication, worsened sharply below 40 runs, suggesting that 50 is close to the minimum size compatible with the Latin square constraints.

Materials and Methods

Materials

The research was conducted at the Munich Center for Computational Sciences, Munich, Bavaria, Germany, between June 2020 and December 2022. Farm-level production records-daily milk yield (kg cow⁻¹), fat percentage, protein percentage, and somatic cell count-were supplied by the British Columbia Milk Marketing Board under a data-sharing agreement. Sixteen farms in the Fraser Valley participated voluntarily; each farm had a herd size between 80 and 200 lactating Holstein cows. Feed-composition data, including batch-level analyses of dry matter, crude protein, neutral detergent fibre, and mineral content, were recorded by farm nutritionists at the start of each milking season (spring and autumn). Weather data for each farm's postal code were obtained from Environment Canada to control for heat-stress effects. The dataset contained 38,412 cow-season records after excluding animals with fewer than 60 days in milk.

Methods

Each of the four feed factors (corn silage proportion, alfalfa hay proportion, canola meal proportion, mineral premix rate) was set at five equally spaced levels spanning the range recommended by dairy nutrition guidelines [10] . The 50-run design was assigned to 16 farms over two seasons using a randomisation schedule whose SHA-256 hash was published in a pre-registration document before the trial started [4] . Milk yield was analysed as a linear mixed model with fixed effects for the four factors and their pairwise interactions, and random effects for farm, season, and farm-by-season. Trimmed means (10% symmetric trim) replaced ordinary means for fat percentage to down-weight outlying observations caused by sampling artefacts. The hazard function for time to production drop-off (defined as the first week a cow's yield fell below 20 kg day⁻¹) was estimated by a Cox proportional-hazards model with feed-factor levels as covariates [3] . Variance components were extracted by restricted maximum likelihood (REML) and expressed as percentages of total variability. Treatment effects dominated the variance structure, accounting for nearly half of total variability. The treatmentby-block interaction was surprisingly large at 18.7%, indicating that farms responded differently to the same feed combinations-a pattern that complicates generalisation across the valley. Higher corn silage proportions were associated with lower hazard of early production drop-off (HR = 0.87 at level 5), while higher canola meal proportions increased the hazard (HR = 1.22). Mineral premix had no meaningful effect on drop-off timing.

Results

Discussion

The convex relaxation produced a 50-run fractional factorial plan with D-efficiency 0.89, which in practical terms means the plan captured 89% of the information that the full 625-run design would have provided. For a multi-farm dairy trial where each run commits an entire herd to a specific diet for months, this reduction from 625 to 50 runs translates directly into saved resources and reduced animal welfare risk [6] .

The Latin square constraints proved important for balancing farm-level confounding. Without them, the Fedorov exchange algorithm tended to concentrate runs on a subset of farms, leaving others under-represented. Imposing equal-replication conditions for each factor within each farm ensured that every farm contributed to every main-effect contrast, a property that strengthened the internal validity of the trial [2] . The treatment-by-block interaction (18.7% of variance) deserves attention. It suggests that optimal feed composition may differ across farms, perhaps because of variation in barn climate, cow genetics, or management style. A follow-up trial could use a split-plot design in which farm is the whole-plot factor and feed composition the sub-plot factor, allowing direct estimation of farm-specific optimal diets [5] . The hazard analysis added a time-to-event dimension that ordinary variance analysis misses. The protective effect of higher corn silage proportion aligns with nutritional theory: corn silage provides glucogenic energy that sustains peak lactation longer [10] . The increase in hazard at high canola meal levels may reflect ruminal protein imbalance, although the confidence intervals were wide enough that this effect should be confirmed in a dedicated feeding trial.

Limitations include the restriction to two milking seasons, which prevents estimation of year effects, and the use of aggregate farm records rather than individual-cow data for some variables. The cryptographic hash step, while useful as a transparency measure, does not guarantee that farms adhered perfectly to assigned diets; spot audits would strengthen compliance monitoring.

Conclusion

This research applied convex optimisation to design a 50-run fractional factorial plan for a four-factor, five-level dairy feeding trial in the Fraser Valley, achieving a D-efficiency of 0.89 under Latin square balance constraints. Variance decomposition showed that feed-factor treatments explained about half of the variability in milk yield, with a noteworthy 18.7% attributed to treatment-by-block interactions. A Cox hazard analysis revealed that higher corn silage proportions delayed production drop-off, while elevated canola meal levels accelerated it.

The sensitivity analysis confirmed that 50 runs sat near the efficiency knee: adding more runs yielded diminishing returns, while reducing the count below 40 broke the Latin square balance. These findings demonstrate that modern optimisation tools can make multi-level factorial experiments feasible in settings where the full factorial is impractical. The cryptographic hashing of the randomisation sequence added a layer of trial integrity that could become standard practice in open agricultural experiments. Future directions include extending the convex framework to split-plot and strip-plot designs, and exploring response-surface designs that allow curvature estimation within each factor.

Acknowledgement

Acknowledgements Funding SourcesThe authors received support from internal research allocations at the Munich Center for Computational Sciences. No external grant number is associated with this work. Institutional SupportThe Department of Mathematical Physics at the Munich Center for Computational Sciences provided access to MATLAB licenses and computing resources. Contributions Not Qualifying for AuthorshipThe authors thank Dr. Hans Becker for helpful comments on the convex relaxation formulation.

References

1. Atkinson AC, Donev AN, Tobias RD. Optimum Experimental Designs, with SAS. Oxford: Oxford University Press; 2007.
2. Bailey RA. Design of Comparative Experiments. Cambridge: Cambridge University Press; 2008.
3. Cox DR. Regression models and life-tables. Journal of the Royal Statistical Society Series B. 1972;34(2):187220.
4. Menezes AA, van Oorschot PC, Vanstone SA. Handbook of Applied Cryptography. Boca Raton: CRC Press; 1996.
5. Fisher RA. The arrangement of field experiments. Journal of the Ministry of Agriculture. 1926;33:503-513.
6. Wu CFJ, Hamada MS. Experiments: Planning, Analysis, and Optimization. New York: John Wiley and Sons; 2009.
7. Box GEP, Hunter JS, Hunter WG. Statistics for Experimenters. New York: John Wiley and Sons; 2005.
8. Dey A, Mukerjee R. Fractional Factorial Plans. New York: John Wiley and Sons; 1999.
9. Kiefer J. Optimum experimental designs. Journal of the Royal Statistical Society Series B. 1959;21(2):272-319.
10. National Research Council. Nutrient Requirements of Dairy Cattle. Washington, DC: National Academies Press; 2001.
11. Fedorov VV. Theory of Optimal Experiments. New York: Academic Press; 1972.
12. Pukelsheim F. Optimal Design of Experiments. Philadelphia: SIAM; 2006.
13. Montgomery DC. Design and Analysis of Experiments. New York: John Wiley and Sons; 2017.
14. Schoen ED, Mee RW. Two-level designs of strength 3 and up to 48 runs. Journal of the Royal Statistical Society Series C. 2012;61(1):163-174.
15. Gilmour SG, Trinca LA. Optimum design of experiments for statistical inference. Applied Statistics. 2012;61(3):345-401.
16. Xu H, Wu CFJ. Generalized minimum aberration for asymmetrical fractional factorial designs. Annals of Statistics. 2001;29(4):1066-1077.
17. Tang B, Deng LY. Minimum G2-aberration for nonregular fractional factorial designs. Annals of Statistics. 1999;27(6):1914-1926.
18. Draper NR, Smith H. Applied Regression Analysis. New York: John Wiley and Sons; 1998.
19. Satterthwaite FE. An approximate distribution of estimates of variance components. Biometrics Bulletin. 1946;2(6):110-114.
20. Kenward MG, Roger JH. Small sample inference for fixed effects from restricted maximum likelihood. Biometrics. 1997;53(3):983-997.
21. Harville DA. Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association. 1977;72(358):320-338.