A Spatial Econometric Analysis of Weather Effects on Milk Production

1. Introduction

The increasing concentration of greenhouse gases (GHGs) in earth’s atmosphere is expected to increase air temperatures, weather variability, and the frequency of weather-related extreme events such as heat waves and freezes [1]. The literature indicates milk production, feed intake, reproductive efficiency and animal health are all sensitive to such changes [2]. Regional effects are expected to be uneven. Understanding impacts of weather condition changes on regional milk production in the U.S. is crucial for the U.S. agricultural sector and attracting much attention from researchers [3,4].

In this study, we use panel econometric models to estimate relationships between weather indicators and milk productivity. Panel data approaches have been widely utilized to analyze climatic effects on productivity as reviewed in Blanc and Schlenker [5]. However, such approaches may suffer from omitted variable biases [6]. In particular, time-varying omitted variables that are common across regions cannot be captured by regional fixed effects. Weather factors, such as wind or solar radiation, vary over time and are spatially correlated. When omitted such variables will lead to spatial correlated error terms. Spatial panel models have been developed to handle such issues [7,8,9,10]. To the best of our knowledge, we are not aware of any study that has applied spatial panel models to examine the extent to which changes in weather condition can affect milk production to reduce omitted variable biases. The rest of the paper is organized as follows. Section 2 reviews related literature. The model structure, assumptions, and data are described in Section 3 and is followed by the model estimation results in Section 4. Conclusions and discussion are provided in Section 5.

2. Literature Review

2.1. Weather and Milk Productivity

Heat stress and the zone of thermal comfort play an important role in livestock production and reproduction [11]. When conditions exceed the upper limit of that zone, livestock exhibit degraded performance [12]. West [13] and Lacetera et al. [14] indicate dairy cows are especially sensitive to heat stress exhibiting reduced feed intake, increased water intake, decreased meat and milk production, and altered birth rates. Ravagnolo, Misztal and Hoogenboom [15] suggest that a one unit increase in heat stress leads to about a 0.2 kg daily milk yield loss compared to a base of 26.3 kg. Mauger et al. [3] estimate climate is causing present-day US Holstein milk production losses of 1.9 percent and project that climate change could cause this to increase to 6.3 percent by 2100 generating a $2.2 billion annual loss.

Heat stress arises from the combination of temperature, air movement, and humidity [11,16,17,18]. The Temperature Humidity Index (THI) is the most commonly used index [19,20] to quantify heat stress. There are several alternative THI formulae as summarized by Bohmanova et al. [16]. The determinants are dew point temperature, wet bulb temperature, relative humidity, wind speed, solar radiation and water vapor [21].

In this study, we use seasonal maximum and minimum THI calculated with the formula of Amundson et al. [22] and apply spatial panel models to relate it to milk production while also including other weather factors. Hammami et al. [23] show milk production exhibits an inverse U-shaped curve with respect to THI with the threshold falling in the range of 62 – 82. In particular, milk productivity is expected to increase at a decreasing rate when the THI increases before reaching a threshold level and then milk yield decreases at an increasing rate as the THI continues increasing.

2.2. Econometric Analysis of Climate Effects on Agricultural Production

There have been a number of panel econometric studies measuring weather impacts on agricultural productivity as reviewed in Blanc and Schlenker [5]. For crop yields, such studies include: Schlenker and Roberts [24], Deschênes and Greenstone [25], Miao, Khanna, and Huang [26], etc. Regionally, Mukherjee, Bravo-Ureta and De Vries [27] examine climatic impacts on southeastern U.S. milk yields, finding that THI has a significant non-linear negative effect. Qi, Bravo-Ureta and Cabrera [28] examine temperature and precipitation effects on Wisconsin dairy farm productivity and find those weather indicators have exerted deleterious effects and would further reduce milk output by 5 to 11% per year between 2020 and 2039. Perez-Mendez, Roibas and Wall [29] estimate effects of weather condition on Spanish milk production considering cow performance and forage supplies. But they find changes in weather mainly affect forage availability and the THI has no significant effect on milk production. Key and Sneeringer [30] examine heat stress in US major milk producing regions, finding it would depress milk production and cause national economic loss.

Although panel models control some forms of omitted variable bias [5,31], they do not control all forms. In particular omitting climate variables that are related across regions (e.g., vapor pressure deficit and wind speed) would bias estimation results [32,33]. This is a concern since climate variables inherently vary across time and are correlated across space [34]. Auffhammer et al. [35] argue that omission of such variables will cause spatially correlated residuals, which will then bias estimation results, and that the situation can be improved by using spatial panel models.

Several studies have employed spatial panel models in studying climate effects on crop yields. Schlenker, Hanemann and Fisher [9] use the Spatial Error Model (SEM), assuming spatially dependent error terms that are uncorrelated with independent variables. Ortiz-Bobea [10] uses the Spatial Durbin Model (SDM) that assumes spatially correlated error terms that are also correlated with independent variables. Herein, this study utilizes both SEM and SDM, analyze their estimation results, and compares model performance with conventional panel data models to explore the feasibility of spatial panel models.

3. Materials and Methods

3.1. Model Specification

A standard panel model [36,37] can be expressed as follows:

$$y=X\beta +\mu$$

(1)

Where in our case $y$ is annual milk production per cow; $X$ is a vector of independent variables (i.e., weather indicators in our study); $\mu$ is an idiosyncratic error that is assumed to be independently identically distributed (iid).

We use two modeling approaches that improve estimations by controlling spatial correlation that arises due to omitted variables. The first is the SEM model of LeSage [38]. SEM assumes the omitted covariates are orthogonal to the independent variables but exhibit spatial dependence. That is in the above model $\mu =\gamma W\mu +\epsilon$ or equivalently,

$${\mu =(I_n-\gamma W)}^{-1}\epsilon$$

(2)

where $\gamma$ is a spatial dependence parameter; ${\mathit{I}}_n$ is an identity matrix;$W$ is a spatial weight matrix that reflects interactions in errors between regions; and $\epsilon$ is an iid error term. In this approach, equation (1) is changed to be:

$$y=X\beta +{(I_n-\gamma W)}^{-1}\epsilon$$

(3)

Second is the SDM model of Cook, Hays and Franzese [39]. Therein the omitted variables are assumed to be correlated with the independent variables. In this case, we assume the error term in (3) is given by

$$\epsilon =X\xi +v$$

(4)

where $\xi$ is the interrelationship coefficient; and $v$ is an iid error term. Then inserting (4) into (3), we obtain:

$$y=\gamma Wy+X(\beta +\xi )+WX(-\gamma \beta )+v$$

(5)

For simplicity of notation, let

$$\tau =\beta +\xi ,\phi =-\gamma \beta$$

(6)

resulting in

$$y=\gamma Wy+X\tau +WX\phi +v$$

(7)

Equation (7) is the expression of the SDM estimating model and includes both the spatially lagged dependent variable and the independent variables. When $\xi =0$ in equation (4) the SDM reduces to the SEM.

Following Belotti, Hughes and Piano Mortari [40], a unifying specification encompassing both cases is:

$$y_{nt}=\alpha +\rho W{y}_{nt}+X_{nt}\beta +W{X}_{nt}\theta +{\varphi }_n+c_t+{\mu }_{nt}$$

(8)

$${\mu }_{nt}=\lambda W{\mu }_{nt}+{\epsilon }_{nt}$$

Where $W$ is a spatial weight matrix; $\rho ,\theta ,\lambda$ are spatial correlation parameters;${\varphi }_n$ and $c_t$ are the individual- and time-specific effects; $u_{nt}$ is a spatially correlated error; ${\epsilon }_{nt}$ is an idiosyncratic error.

In (8) if $\rho$ and $\theta$ are all equal to zero, the model reduces to the SEM case with the spatial interaction in the error term as denoted by $\lambda W{\mu }_{nt}$, which means that milk production in one region could be affected by unobserved factors in neighboring regions.

On the other hand, if $\lambda =0$, the model reduces to the SDM case, and $\beta$ and $\theta$ are coefficients expressing the direct and indirect (spatial spillover) effects of the independent variables. Finally, if $\rho =\lambda =\theta =0$, this reduces to a common panel model without spatial interaction.

In the spatial models, the matrix $W$ reflects regional proximity. In this case, we follow Fischer and Getis [41] and set the elements to $1/n$ for each of n neighboring states and zero for non-neighboring states. We also follow Taha et al. [42] and use the queen criterion to identify neighboring states.

3.2. Data

The above non-spatial and spatial panel models are estimated with the dependent variable being annual milk production per cow (MPC), and the explanatory variables being state level seasonal values of the palmer drought severity index (PDSI), total precipitation, and maximum and minimum THI. We use a dataset from 1950 to 2022 for the 48 contiguous US states. This yields 3504 observations and a balanced panel. Statistics computed over the data are shown in Table 1.

Data sources and manipulations are as follows:

Historical annual milk production: data on annual milk production per cow are drawn from United States Department of Agriculture’s National Agricultural Statistics Service Quick Stats database¹, and are measured in pounds per cow in a year (lbs/head).

Historical climatic data: state level daily palmer drought severity index (PDSI), total precipitation, minimum and maximum temperature, and minimum and maximum relative humidity (RH) data are drawn from National Oceanic and Atmospheric Administration (NOAA) [43]. These data are aggregated into 4 seasons (spring-March to May, summer-June to August, fall-September to November, and winter-December to February). Over these periods, PDSI is the average value, and the more negative the value the drier the climate is with positive values being wetter; precipitation data are the three-month total. Daily minimum temperature and maximum RH are used to calculate daily minimum THI (Min THI), whereas daily maximum THI (Max THI) are calculated using maximum temperature and minimum RH [44,45], following the THI formula of Amundson et al. [22] and Mader, Davis, and Brown-Brandl [17]. Seasonal Min THI and Max THI are computed by averaging the values across the days in the three-month period. Temperatures are measured in degrees Celsius, and precipitation is measured in inches. We use the data from 1950-2014 to execute model estimation and reserve data from 2015-2022 for out-of-sample testing.

In our estimation, quadratic terms for the weather variables are included, following previous studies [46,47,48], as milk yield tends to be non-linearly correlated with a given weather variable. A linear and a quadratic time trend are also added as independent variables to capture "the effect of technological progress". The dependent variable – MPC is transformed into log form (as recommended by Mauldon [49] and Mosheim [50]), which helps stabilize the variance in the data [51].

4. Results

Estimations are done using conventional and spatial panel regression models. We estimate a pooled OLS, a fixed effects panel, spatial panel models (i.e., SEM and SDM).

4.1. Results without Spatial Interactions

Pooled OLS and fixed effects panel model results are summarized in Table 2. Therein most weather variables exhibit a U-shaped effect with a significant quadratic term. The fixed effects model fits better as the null hypothesis of no state-specific fixed effects is rejected at the 1% level using the standard F test. The climate effects in the pooled OLS model are generally larger, indicating neglecting regional fixed effects overstates the climate impact. The smaller results under the panel model likely occur since its fixed effects terms capture important regional omitted variables such as soils, elevation, topography, etc.

In the fixed effects panel model, the summer max THI exhibits an inverse U-shaped effect, meaning that holding other factors constant, MPC will increase as the THI value increases up to a threshold value, and then if the summer max THI increases beyond that threshold level MPC will decrease. This is generally consistent with findings from previous research [52].

4.2. Spatial Panel Model Results

The spatial panel models are relevant when there is spatial dependence in the error terms. We test this using Pesaran’s CD test [53,54] computed over the fixed effects panel model residuals. As a result, the no cross-regional dependence null hypothesis is rejected at the 1% level. Hence, it is appropriate to move on to estimation results from the spatial panel models.

The SEM and SDM estimation results are shown in Table 3. For the SEM, the spatial parameter λ is positive (0.67) and significant at the 1% level, implying a strong spatial dependence in the unobserved factors. The effects of a random shock in a specific region that transmit to its neighborhood can be estimated by multiplying λ with the corresponding spatial weights.

Now turning to the SDM results, the spatial autoregressive parameter ρ (0.62), which reflects the average intensity of the spatial intercorrelation [38], is also significant and indicating a strong spatial dependence. The spatially lagged independent variables ($WX$) measure the indirect impacts on MPC in one region arising from changes in variables $X$ from its neighbored regions [55]. Following LeSage and Dominguez [56] and Sarrias [57], weather influences can be split up considering the direct in-region effects of $X$ and the indirect effects translated in from adjacent regions. In particular, equation (7) can be rewritten to equation (9):

$$y_{nt}={\sum }_{k=1}^K{\left(I_n-\gamma W\right)}^{-1}\left(I_n{\tau }_k+W{\phi }_k\right)x_{k,nt}+{\left(I_n-\gamma W\right)}^{-1}{\nu }_{nt}$$

(9)

and the effects of the independent variables can be computed as

$${\mathit{S}}_k=\frac{\partial y_{nt}}{\partial x_{k,nt}}={\left({\mathit{I}}_n-\gamma \mathit{W}\right)}^{-1}\left({\mathit{I}}_n{\tau }_k+\mathit{W}{\phi }_k\right)$$

$$={\left({\mathit{I}}_n-\gamma \mathit{W}\right)}^{-1}\left[\begin{array}{ccc}\begin{array}{cc}{\tau }_k& w_{12}{\phi }_k\\ w_{21}{\phi }_k& {\tau }_k\end{array}& \cdots & \begin{array}{c}{w}_{1n}{\phi }_k\\ w_{2n}{\phi }_k\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{cc}{w}_{n1}{\phi }_k& w_{n2}{\phi }_k\end{array}& \cdots & {\tau }_k\end{array}\right]$$

Where $n$ identifies region, $k=1,\dots ,K$ denotes the set of independent variables, and ${\mathit{S}}_k$the effect on production across regions of variations in the $k^{th}$ independent variable. In turn, the average direct effect is the average of the diagonal elements of ${\mathit{S}}_k$, and the average total effect is equal to $n^{-1}{\iota }_n^{\text{'}}{\mathit{S}}_k{\iota }_n$, which ${\iota }_n$ is a vector of 1’s. The indirect effect is the total effect minus the direct effect.

For our estimation the empirical SDM direct and indirect impacts are shown in Table 4, along with the SEM results for comparison. Overall, the SDM estimates show that more weather variables exert significant effects, relative to the SEM results. Summer max THI again exhibits an inverse U-shape total effect on MPC; thereby milk yield increases as the summer max THI increases toward a threshold value, and then yield declines as this THI is higher than the threshold moving extreme heat stress, which coincides with the findings in Hammami et al. [23]. We also find a significant indirect effect of summer max THI, which represents the spatial spillover impacts on MPC in a state from changes in summer max THI of its nearby states [10].

Additionally, the linear and quadratic time trend variables were included here to capture the effects of factors like technological improvements changing over time, such as developments in animal genetics, nutrition and animal management. In both spatial models and the pooled OLS and non-spatial fixed effects panel model, the linear time trend is significant and positive showing increasing MPC over time, while the quadratic term is negative, showing diminishing technological progress over time.

These results can be illustrated graphically. Figure 1 displays the relationship between a given weather factor with the MPC when holding other variables constant. Figure 1a shows the percentage changes in MPC as the summer max THI increases by one unit based on the historical minimum (64.8). The summer max THI has an inverse U-shaped effect with a threshold value of 72. Above that threshold, the increase in the summer max THI will decrease MPC. However, if the THI is below this threshold value, its increase will positively affect the MPC, because in that lower THI range, warmer environment would benefit livestock production as found in Du Preez et al. [58] and Correa-Calderon et al. [59]. Nationally, most regions experienced MPC reduction induced by heat stress as the comparison between the historical summer max THI observations and the threshold shows more than 2/3 of past observations exceed the threshold value of 72.

The effects of other weather factors are also plotted in Figure 1. Fall precipitation also exhibits an inverse U-shaped relationship with MPC (Figure 1b), that is, MPC increases at a decreasing rate as fall precipitation rises up to the threshold value of 14 inches. Beyond that threshold value, MPC declines at an increasing rate. Generally, since about 90% of the fall precipitation observations are smaller than the threshold, to some extent, more precipitation would increase MPC perhaps reflecting added forage availability [60,61]. For winter PDSI (Figure 1c), its total effect also shows an inverse U-shape with a threshold value of 0. This implies as PDSI rises indicating less drought, MPC tends to increase, but as PDSI continues increasing to be higher than the threshold value of 0, then MPC falls as environmental condition gets wetter. These results agree with previous literature [62,63].

4.3. Model Comparison

Now we compare model performance by examining the in and out of sample model fit using the root mean squared error (RMSE) as the criterion [64]. The in-sample RMSE is computed with the data from 1950 to 2014, assessing how effective the models are in reproducing data, whereas the out-of-sample RMSE, calculated by the holdout dataset from 2015 to 2022, is to evaluate the forecasting accuracy of the models on a new dataset of independent variables [65]. The results in Table 5 show that the estimates from the spatial models outperform that from the conventional panel data models.

Next, we statistically test differentiating parameters using a Wald test. First, we test whether the SEM spatial dependence (λ) factor is significantly different from zero and find it is at the 1% confidence level. In turn, we conclude the SEM model is more appropriate than the fixed effects panel model and that it is important to incorporate spatial correlation.

Next, we test whether $\phi +\gamma \tau =0$, which if so, means that SDM reduces to SEM [66]. This is also rejected at the 1% level, which indicates that SDM performs better than SEM. Thus, it is important to incorporate the correlation between the independent and the omitted variables.

5. Conclusions and Discussion

Our estimation results show changes in weather conditions influence the U.S. milk production (on a per cow basis). Summer heat stress, as measured by summer max THI, exerts a significant effect on milk production. They are estimated to have an inverse U-shaped relationship with a threshold value of 72, when holding other variables constant; below this threshold, milk yield increases as the summer max THI rises, but when this THI continues increasing to be higher than that threshold value, milk productivity declines, which are aligned with findings by Renaudeau et al. [52]. Fall precipitation and winter drought conditions also exert an inverse U-shaped effect with threshold values of 14 inches and 0 on the PDSI scale, respectively. In essence, regarding such inverse U-shaped relationship, an extreme weather condition, such as moving toward too much or too little rain, will negatively affect milk production.

On the methodological side, we use several model forms in an effort to reduce omitted variable bias and improve model fit. The results strongly support spatial panel models, particularly the SDM, that can control spatial correlation that arises due to omitted weather variables. It should be noted that the spatial panel models are designed to reduce the effect of omitted spatially pervasive forces but does not eliminate that bias. Therefore, in subsequent research, it would be good to explore improved model specifications and include more climate related variables.