lactationcurve.characteristics.centering_date_method

  1"""
  2Centering date method for accumulated lactation yield.
  3
  4Purpose
  5-------
  6This module implements the centering date method to compute cumulative milk yield
  7from test-day records.
  8This method is also known as the McDowell mehod for centering the testing day.
  9This approach was the standard method used by U.S. DHIA before 1969.
 10
 11Method Summary
 12--------------
 13The centering date method (CDM) estimates lactation yield by assuming that each
 14test-day milk yield is representative for the calender month it was taken in.
 15It estimates total milk yield over a lactation by taking a milk yield measurement
 16each month and multiplying it by the number of days in that month.
 17If a cow calved before the start of the month,
 18the extra days are multiplied by the first test-day milk yield. If a cow goes
 19dry after the end of the month, the extra days are multiplied by the last
 20test-day milk yield.
 21So each testing interval has the same amount of days as days in that month,
 22whereby the testing day is considered as the middle of the testing period.
 23The total lactation yield is then calculated by multiplying each test-day yield
 24by the length of its corresponding centered interval and summing these
 25contributions across the lactation.
 26
 27
 28Formula:
 29
 30MY_total = (d₁ + ℓ₁/2) · y₁ + Σ ℓᵢ · yᵢ (i = 2 … n−1) + (L − dₙ + ℓₙ/2) · yₙ
 31
 32Where:
 33- ``d``: days in milk of the test day
 34- ``i``: index of the test day in the lactation
 35- ``yᵢ``: milk yield measured in the 24 hours of the test day i
 36- ``ℓ``: number of days in the month that the measurement was taken
 37- ``L``: lactation length (e.g., 305 days or the days in milk that the cow went dry)
 38- ``n``: number of test days in the lactation
 39- ``MY_total``: total cumulative milk yield over the lactation
 40
 41
 42
 43Key Entry Points
 44----------------
 45- ``centering_date_method``: Computes cumulative lactation milk yield for
 46    a group of lactations, each identified by a unique ``TestId``
 47    and returns the results as a DataFrame with a ``TestId`` and ``LactationMilkYield`` column.
 48- ``centering_date_method_single_lactation``: Computes cumulative lactation
 49    milk yield for a single lactation.
 50
 51
 52Column Flexibility
 53------------------
 54The function accepts several case-insensitive column name aliases and can
 55create a default ``TestId`` if one is missing. Recognized aliases:
 56
 57- Days in Milk: `["daysinmilk", "dim", "testday"]`
 58- Milk Yield: `["milkingyield", "testdaymilkyield", "milkyield", "yield"]`
 59- Test Id: `["animalid", "testid", "id"]`
 60
 61It is also possible to provide your own column names so the function
 62can be applied to dataframes with different column naming conventions.
 63
 64The main function returns a DataFrame with columns: ``["TestId", "LactationMilkYield"]``.
 65
 66Notes
 67-----
 68- Units: DIM (days in milk) is measured in days,
 69    and milk yield can be in kg or lb. The
 70    output stays in the same unit as the input.
 71- Records with ``DIM > lactation_length`` are excluded before computation.
 72- This method's main strength is its ease of use and simplicity.
 73- Its main disadvantages are that it does not account for the shape of the
 74    lactation curve, which can lead to overestimation or underestimation
 75    of total yield,
 76    especially for lactations with few test days or irregular patterns.
 77    Outliers in test-day records can also have a large influence on the final
 78    result.
 79- This deviation is most prevalent at the beginning and end of the lactation
 80    as it ignores the increase in milk yield at the start of lactation
 81    and the decrease at the end of lactation.
 82    Everett et al. (1968) has proposed correction factors to account for this,
 83    but they are not implemented in this method.
 84    The ICAR guideline allows for adjustment of the first recorded value when
 85    computing the accumulated record.
 86    If the first recording occurs within 14 days from calving,
 87    then no adjustment is required to the first recorded value when computing
 88    the accumulated record. If the first recording occurs 15 to 95 days from
 89    calving, then an adjustment procedure may be applied.
 90- Months are of inequal length so the cummulative milk yield is biased
 91    by when the cow freshens, e.g. a cow that has her peak yield in January
 92    will likely have a higher cumulative milk yield.
 93    according to this method than a cow that has her peak yield in February,
 94    even if they have the same total milk yield over the lactation.
 95- Because the method assumes a once per month testing scheme,
 96    you cannot use it with a testing scheme that has multiple tests per month.
 97- This method works best with a very consistent once per month testing scheme.
 98    For instance always take the first day of the month so that the interval
 99    length is always the same.
100- If a month is missing, the milk yield for that month is estimated as the average
101    of the nearest test before and the nearest test after.
102- If the first test day is missing, the second test day is used to estimate
103the milk yield for the first month.
104- This method can also be applied to bimonthly or trimonthly testing schemes,
105    but the accuracy rapidly decreases as the number of test days per lactation
106    decreases.
107    This has been shown by Everett et al. (1968) reporting errors up to 28%
108    and is also true for the test interval method.
109-The code takes into account leap years if dates include years.
110    The milk yield of February will than be multiplied by 29 instead of 28.
111- This method is different from the test interval method **only**
112because it uses the number of days in the month as the interval length
113instead of the actual number of days between test days.
114Otherwise they mathematically yield the same result.
115
116
117
118Reference
119---------
120U.S. Department of Agriculture, Bureau of Dairying. (1925).
121Cow Testing Association Letter No. 1. Washington, DC.
122
123R.W. Everett, B.T. McDaniel, H.W. Carter,
124Accuracy of Monthly, Bimonthly, and Trimonthly Dairy Herd Improvement Association Records,
125Journal of Dairy Science,
126Volume 51, Issue 7,
1271968,
128Pages 1051-1058,
129ISSN 0022-0302,
130https://doi.org/10.3168/jds.S0022-0302(68)87123-1.
131
132
133---------
134Author: Meike van Leerdam
135Date: 07-31-2025
136Last update: 22-May-2026
137"""
138
139from typing import cast
140
141import numpy as np
142import pandas as pd
143
144from lactationcurve.preprocessing import standardize_lactation_columns
145
146
147def centering_date_method_single_lactation(
148    lactation_df: pd.DataFrame,
149    lactation_length: int = 305,
150    lactation_start_grace_period: int = 0,
151) -> float:
152    """Compute total lactation milk yield using the centering date method.
153
154    Args:
155        lactation_df (pd.DataFrame): Input DataFrame that contains all milk
156            recordings of a single lactation
157            with at least TestCalendarDate,DaysInMilk and MilkingYield columns.
158            Column names can be provided explicitly or matched via known aliases.
159        lactation_length (int): Lactation length used to calculate cumulative
160            production. The default is 305 days.
161            Records with DIM > lactation_length are excluded.
162            If dry dates are known it is recommended to use the actual lactation
163            length instead of the default 305 days.
164        lactation_start_grace_period (int): Number of days after calving that
165            are not considered for cumulative yield becuase the colostrum is
166            not representative for the lactation. The default is 0 days.
167
168    Returns:
169        float: The computed total milk yield over the lactation for a single cow.
170
171    """
172    df = lactation_df.sort_values("TestCalendarDate")
173
174    # One yield per calendar month (if a month has several tests, take their mean)
175    month = df["TestCalendarDate"].dt.to_period("M")
176    y = cast(pd.Series, df.groupby(month)["MilkingYield"].mean())
177
178    # Full monthly range from first to last test; missing months become NaN
179    full_range = pd.period_range(y.index.min(), y.index.max(), freq="M")
180    y = y.reindex(full_range)
181
182    # Missing month = average of the nearest test before and the nearest test after
183    y = y.fillna(y.ffill().add(y.bfill()).div(2))
184
185    # Weights: days in each month, with the first and last month adjusted
186    w = np.array([period.days_in_month for period in full_range], dtype=float)
187    d_first = df["DaysInMilk"].iloc[0] - lactation_start_grace_period
188    d_last = df["DaysInMilk"].iloc[-1]
189    w[0] = d_first + w[0] / 2
190    w[-1] = lactation_length - d_last + w[-1] / 2
191
192    total_milk_yield = (w * y).sum()
193    return total_milk_yield
194
195
196def centering_date_method(
197    df: pd.DataFrame,
198    days_in_milk_col: str | None = None,
199    milking_yield_col: str | None = None,
200    test_id_col: str | None = None,
201    test_calender_date_col: str | None = None,
202    default_test_id: int = 0,
203    lactation_length: int = 305,
204    lactation_start_grace_period: int = 0,
205) -> pd.DataFrame:
206    """Compute total lactation milk yield using the centering date method for each test id.
207
208    The method applies:
209
210
211    Args:
212        df (pd.DataFrame): Input DataFrame with at least DaysInMilk and
213            MilkingYield columns, plus an optional TestId column. Column names
214            can be provided explicitly or matched via known aliases.
215        days_in_milk_col (str | None): Optional column name override for
216            DaysInMilk.
217        milking_yield_col (str | None): Optional column name override for
218            MilkingYield.
219        test_id_col (str | None): Optional column name override for TestId.
220        test_calender_date_col (str | None): Optional column name override for
221        TestCalendarDate.
222        default_test_id (int): Value used to create a default TestId column if
223            one is missing.
224        lactation_length (int): Lactation length used to calculate cumulative
225            production. The default is 305 days.
226            Records with DIM > lactation_length are excluded.
227        lactation_start_grace_period (int): Number of days after calving that
228            are not considered for cumulative yield becuase the colostrum is
229            not representative for the lactation. The default is 0 days.
230
231    Returns:
232        pd.DataFrame: Two-column DataFrame with
233            - "TestId": identifier per lactation,
234            - "LactationMilkYield": computed total milk yield over the
235              specified window.
236
237    Raises:
238        ValueError: If required columns (DaysInMilk or MilkingYield) cannot be found.
239        ValueError: If there are multiple measurements in the same calendar month.
240        ValueError: If there are fewer than two data points per TestId.
241
242    """
243
244    # Standardize columns and filter DIM <= lactation_length
245    df = standardize_lactation_columns(
246        df,
247        days_in_milk_col=days_in_milk_col,
248        milking_yield_col=milking_yield_col,
249        test_id_col=test_id_col,
250        test_calender_date_col=test_calender_date_col,
251        default_test_id=default_test_id,
252        max_dim=lactation_length,
253        lactation_start_grace_period=lactation_start_grace_period,
254    )
255
256    result = []
257
258    # Iterate over each lactation
259    for lactation in df["TestId"].unique():
260        lactation_df = pd.DataFrame(df[df["TestId"] == lactation])
261
262        if len(lactation_df) < 2:
263            print(f"Skipping TestId {lactation}: not enough data points for interpolation.")
264            continue
265
266        # Check that the date column is datetime and convert it if needed.
267        if not pd.api.types.is_datetime64_any_dtype(lactation_df["TestCalendarDate"]):
268            lactation_df["TestCalendarDate"] = pd.to_datetime(lactation_df["TestCalendarDate"])
269
270        # Check that there is no more than one test per calendar month
271        month_counts = lactation_df["TestCalendarDate"].dt.to_period("M").value_counts()
272
273        if (month_counts > 1).any():
274            duplicated_months = month_counts[month_counts > 1].index.tolist()
275
276            raise ValueError(
277                f"TestId {lactation}: multiple measurements found in "
278                f"calendar month(s) {duplicated_months}. "
279                "The McDowell centering date method requires maximum one test per calendar month."
280            )
281
282        # Calculate total milk yield for this lactation
283        total_milk_yield = centering_date_method_single_lactation(
284            lactation_df,
285            lactation_start_grace_period=lactation_start_grace_period,
286            lactation_length=lactation_length,
287        )
288
289        result.append({"TestId": lactation, "LactationMilkYield": total_milk_yield})
290
291    return pd.DataFrame(result, columns=pd.Index(["TestId", "LactationMilkYield"]))