Abstract

Statistical modeling is fundamental in understanding and predicting the variability of climate change and its effect on cereal crop yield in Northern Ghana. Most often, there is little or no data to study these relationships within the Northern Ghana enclave. The Dynamic Bayesian Network Modeling was employed to model the relationship using data on some selected cereal crops (Maize, Rice, and Sorghum) yield from Ghana’s Ministry of Food and Agriculture, and some climate variables (rainfall, sunshine, temperature, humidity, and windspeed) from Ghana Meteorological Agency, for Northern Ghana for a thirty-one-year period. The results showed that climate variables are significantly related to yield at the intra-slide and inter-slide levels and that there exist positive and significant (link strength values > 1) relationships between cereal crop yield and climate change variables. The results imply that, with data on current climate variables under the study, the next cereal crop yield can be predicted with reduced uncertainty and better forecasting ability while holding other factors constant. On the basis of the analysis, the study affirmed the use of the Dynamic Bayesian Networks and highly recommended the extension of the study to include other variables such as vegetative cover, cultivated area, and chemical use, when data becomes available.

1. Introduction

It is observed that the six ecological zones in Ghana that provide livelihood support (forestry and agricultural) and other critical natural resources for rural communities and the country’s economy are not spared the clutches of climate change ¹. Agriculture, being a major source of employment in Ghana suffers the most from climate change; and predominantly affected by the distribution of rainfall. The increasing rate at which rainfall patterns vary result in increase in the risk linked to farming, since it makes prediction very difficult ^{2, 3}.

The Ghanaian economy to some extent, relies on sectors that are sensitive to climatic conditions, mainly agriculture, energy and forestry. The inadequate use of irrigation resources and high dependence on unfavorable climatic conditions for good harvest tend to introduce significant volatility in the standard of living of people. Climate change poses a serious socioeconomic risk to Ghana, especially in the Northern sector which has a long dry season of about seven months and a shorter rainy season of about five months with mostly intermittent drought and floods ⁴. The Northern sector is also observed to be the most agricultural part of Ghana with a high percentage of subsistent food crop farmers, and manifestations of climatic instability in the form of floods, late rains before planting season, and persistent droughts during the planting season ⁵.

Northern Ghana also contributes immensely to the production of rice and maize in the country and is a primary source of sorghum and millet in Ghana ^{6, 7, 8}. Cereals such as maize, rice, sorghum, and millet are sources of major staples in Northern Ghana in particular and the country as a whole ⁶, besides its other commercial uses in areas such as the poultry and brewery industries ⁹. The decision to conduct the study in Northern Ghana is therefore worthwhile.

Further, the effects of climate change on food crop production is getting worse ^{4, 10}, and it is even affecting the achievement of the United Nation’s Sustainable Development Goal Two; which seeks to end all forms of hunger, achieve food security, improve nutrition, and promote sustainable agriculture by 2030 ^{11, 12}. This makes every contribution towards enhancing the efforts at ensuring the sustenance of food production to feed the rising population and the economy, especially in the Northern sector of the country, worthwhile.

The socioeconomic importance of the agricultural sector and cereal production for that matter, in the country cannot be overemphasized, and for that matter, steps towards mitigating the impact of climate change in the Northern sector is crucial. To be able to effectively implement any mitigation scheme in this respect, climate change variability and its effects, especially on cereal crop yield, need to be understood and be able to be predicted. This makes it imperative to include statistical modeling in the process. The question on what the relationship between cereal crop yield and climate change in Northen Ghana arises. The attempt at finding an answer to this question results in the quest to investigate and model the impact of climate change on the yield of some selected cereal crops in Northern Ghana. This is to enable easy and accurate assessment of climate change impacts on cereal crop yield, to enhance understanding of the situation ahead of time, and to inform mitigation actions. It will further contribute to literature on climate change and its impact on cereal crop yield in the area, and will also show the appropriateness of the application of the Dynamic Bayesian Networks in the study.

Climate Change here is, a change of climate that is attributed directly or indirectly to human activity that alters the composition of the global atmosphere and which is in addition to natural climate variability observed over comparable periods ¹³.

Some climate-related research works have been done in the study area. Dumenu and Obeng (2016), observed that the most experienced climate change impacts in the most vulnerable ecological zones (Sudan and Guinea Savanna zones) are prolonged drought, erratic rainfall, reduction in crop yield and frequent flooding ¹⁴. Antwi-Agyei et al. (2012) found the Regions in Northern Ghana, as the most vulnerable to drought. These areas are also observed to have high temperatures and low levels of rainfall ¹⁵. Prolonged drought and erratic rainfall are particularly seriously affecting the rural dwellers in these areas, who rely solely on rains for farming as their occupation ³. Climate change is a danger to agriculture, food security, and national development, and as a result, there is a need to plan adaptation and mitigation actions that will help reduce its effects.

Studies have also been conducted to assess climate change impacts, vulnerability, and adaptation strategies in Ghana ^{16, 17, 18}; impacts on livelihoods and adaptation strategies of rural communities ¹⁹. There are however, little or no sufficient studies on the relationship between climate change and cereal crop yield in Northern Ghana. Further, there is also the issue of inadequate data on cereal crop yield from the study area, hence making the application of a Bayesian approach most suitable.

Bayesian network (BN) models are a kind of probabilistic graphical models that have become very popular among practitioners mainly due to the powerful probability theory involved, which makes them able to deal with a wide range of problems ²⁰. In their work on groundwater quality assessment using data clustering, ²¹ adopted the mixture of truncated exponentials (MTEs) framework in the design and application of a probabilistic clustering, based on hybrid Bayesian networks ²². ²⁰ also applied the Bayesian Networks in their work on climate change assessment using the evaluation of its impact on the olive system in Andalusia. Their results predicted that climate change will lead to changes in the territorial distribution of olive crops. They also observed from their results that Bayesian networks are a powerful tool that allows dealing with both discrete and continuous data ²³.

Dynamic Bayesian Networks (DBNs) are extensions of Bayesian networks to model dynamic processes and consist of a series of time intervals that present the states of all variables at a given time and thus represent the evolution of a process over time ^{24, 25}. They can be seen as a generalization of Markov Chains and Hidden Markov Models since they represent a space of states in a factorized way instead of as a single discrete random variable and can be classified as directed and probabilistic ²⁶. DBNs have been applied in some areas such as pattern detection of human brain behavior ²⁷, speech recognition ^{28, 29}, medical diagnosis ^{30, 31, 32}, visual activities and crime risk analysis ^{30, 33, 34}, sensor validation, client analysis and video object tracking ^{30, 34} to capture their dynamic behavior in the modeling process.

2. Materials and Methods

The study covers climate change variables and cereal crop yield in Northern Ghana, which comprises five administrative regions, namely, the Northern, Savanna, North East, Upper East, and Upper West regions. It lies between latitudes 8°N and 11°N and has a land area of 97702 km², with rainfall characterized by a long dry period of about seven months from October/November to April/May with no appreciable rainfall level, and a rainy season from May to October with adequate levels of rainfall. Annual precipitation ranges between 400 and 1200mm. Agriculture is about 95% purely rain-fed and employs about 70% of the population. The region is one of the most agricultural parts of Ghana with a high percentage of subsistent food crop farmers; covering about 36% of the nation’s cereal (maize, rice, and sorghum) yield over the last three decades ³⁵. The Region is sensitive and vulnerable to climate change and unpredictable climatic variability, especially erratic rainfall ³⁶.

The study relied on secondary data from state institutions mandated to collect and handle such data. Data on climate variables (Rainfall, Temperature, Hours of Sunshine, Relative humidity, Wind speed) for Northern Ghana, was taken from seven weather stations in the North (Bole, Salaga, Tamale, Wa, Yendi, Walewale, and Bawku) from 1992 to 2022 was sourced from the Ghana Meteorological Agency (GMA) through its Client Service Office. Data on Annual Cereal Production figures was also sourced from the Statistics, Research and Information Directorate (SRID) of Ghana’s Ministry of Food and Agriculture (MoFA) for the period 1992 to 2022 for the study area. Data was standardized to help eliminate the effects of other influencing factors and also ensure standardized results for comparability and easy interpretability. Besides, since climate data are measured in different units of measurement, the scaled data is appropriate to synchronize for analysis.

Bayesian network (BN) is a statistical multivariate model for a set of random variables which is defined in terms of two components, namely, the qualitative component, comprising a directed acyclic graph (DAG) of random variables (vertices,), linked by edges between them representing the existence of statistical dependence relationships between these variables in the model; and quantitative component, which is the conditional probability distribution for each variable , given its parents in the graph expressed in conditional probability tables (CPTs) (for discrete variables) or probability functions (for continuous variables) ^{20, 37}. The qualitative component allows Bayesian Network models to be easily understood by experts in other fields who are not familiar with the model's mathematical context; and also allows that, with no mathematical calculation involved, the variable(s) that are or are not relevant for a certain one or more variables can be known, and help to simplify the joint probability distribution of the variables necessary to specify the model ^{38, 39}. Bayesian Networks therefore provide a compact representation of the joint probability distribution over all the variables, defined as the product of the conditional distributions attached to each node, so that

(1)

where is a set of parents of variable , according to the structure of the directed acyclic graph ³⁷.

Various Bayesian Networks have been developed and used in modeling in different situations and on various data sets. Some of these include the Naïve Bayes (NB), Tree Augmented Naïve Bayes (TAN), Mixture of Truncated Exponentials (MTE) for regression ^{40, 41, 42, 43, 44, 45} and Additive Bayesian Network (ABN) ^{46, 47, 48}.

Dynamic Bayesian Networks (DBN) is a type of probabilistic graphical model that extends the static Bayesian Network structure to instances of time series, in which time series data are used for formulating causal relationships among random variables. A DBN is also comprised of a directed acyclic graph defining its structure and a set of parameters that define the probabilistic associations of the variables. In the case of DBNs, time is discretized into time slices that represent consecutive instances.

If a set of variables in a single time slice , and for some horizon , we can define the joint probability distribution of the network as

(2)

where represents the probability distribution of a set of nodes , and represents a set of parent notes of in the graph. Here, is the most recent instance.

In this setting, it is possible that nodes will have parents in previous time slices and for that matter variables in all instants have to be taken into consideration in computing the joint probability distribution as in Eq. (2). In this case, a common assumption is to suppose the DBN to be first-order Markovian ^{49, 50}. This assumes that the future state of the system is independent of the past given the present, that is, .

The Dynamic Gaussian Bayesian Networks (DGBN) is an extension of the Gaussian Bayesian Networks (GBN) model to a dynamic case, where the temporal component of a Dynamic Bayesian Network (DBN) is considered. In this context, time is discretized into time slices in a given period. There is a static GBN in each time slice which has the usual intra-slice parents, and at the same time has parents in the previous slices, and is also connected to the next slice. All the time slices are accounted for by the joint probability distribution until a certain horizon :

(3)

where is the vector of all the nodes in a time slice and . Eq. (3) requires that all the previous time slices should be considered in computing the product. The Markov assumption ⁵⁰ is used to simplify this situation.

Definition 2.5.1. A Dynamic Bayesian Network fulfills the Markov assumption if, for all , the future time is independent of the past given the last time slices. This is known as the Markovian order of the network ^{49, 50}.

The Markovian order determines the number of time slices required to assume with certainty that the present is independent of the past. The resulting joint probability distribution obtained when a Markovian order one is applied to the network will be

(4)

When the temporal component is combined with the GBN structure, DGBNs are obtained, specified by linear Gaussian CPDs in each time slice. To calculate the joint density over time the nodes in the current time slice have to be taken into consideration together with the parent nodes in the previous time slice only.

It is a common practice to assume that the framework of the network in a DBN is homogeneous over time. This implies that the structure remains constant and not independent on . It means the inter-slice and intra-slice arcs and the parameters are reproduced whenever the network is unrolled. If a Markovian order one DGBN is created, it holds that given the present time slice, the future is independent from the past, and the whole Dynamic GBN can be represented by the first two time slices. This implies also that only inter-slice arcs are allowed from to . As the Markovian order increases, more arcs can emerge from earlier lags to the present. If it is detected that the state of the system depends on more than one temporal lag, the Markovian order of the network can be increased at the expense of greater complexity in learning its structure, its parameters, and for performing inference ^{49, 51}.

The structure of the DGBN is learned to appreciate how the present is affected by the past values of the variables by adapting the dataset to the time discretization of DBNs and using the appropriate learning algorithm. Three structure learning algorithms are identified in learning GDBN. These are the dynamic max-min hill-climbing (DMMHC) algorithm ⁵², the binary particle swarm optimization (PSO) algorithm ⁵³, and the natural number order invariant encoding PSO algorithm ⁴⁹.

The structural learning in this study was done with our data set and the dynamic max–min hill climbing (DMMHC) algorithm ⁵⁴,. The DMMHC is an extension of the max-min hill climbing (MMHC) algorithm ⁵⁴, a hybrid method that searches the space of possible structures with a local search and then directs the arcs and scores the resulting networks with an evaluation criterion. The algorithm begins by performing a local search of the structure around the possible parents and children of each node, by determining the Markov blanket of the node. The group of nodes that make a node conditionally independent from the other nodes is the Markov blanket of that node. There are many estimates of conditional independence with continuous variables and given data. One of such is the exact test for Pearson’s correlation coefficient:

(5)

where is the partial correlation coefficient of and given the set of variables and represents the number of nodes in the set .

If is the set of all nodes in the graph, and , then what is essential is to find the set of nodes that makes a node independent from all nodes . When the test is not significant, an undirected arc between 𝑋 and all the nodes in 𝐙 is added to help find the local structure of the nodes. The skeleton of the network is obtained by combining the local structure of all the nodes found, which constrains the step of scoring and finding the best structure ⁵⁵. That is, to direct the arcs in the network, only edges that are present in the skeleton of the net are allowed to be included, and the resulting structures are scored with the Bayesian information criterion (BIC) ⁵⁶. Once the local structure of all the nodes is found, we obtain the skeleton of the network by combining them, which constrains the step of scoring and finding the best structure ⁵⁷. Thus, to direct the arcs in the network, only edges that are present in the skeleton of the net are allowed to be included, and the resulting structures are scored with the Bayesian information criterion (BIC) ⁵⁸; which measures how well a model fits the data, and at the same time, penalizes complexity as measured by the number of parameters in the model. The model parameters are then fitted using maximum likelihood estimation, upon determining the structure of the network.

The exact inference method is implemented in this study by computing conditional probabilities in the equivalent multivariate Gaussian. It consists of converting the DGBN to its equivalent multivariate joint Gaussian , where are the objective nodes and are the instantiated nodes. Once the evidence is instantiated, the objective node is marginalized. When the network is translated to a mean vector and a covariance matrix , the posterior mean and variance of the objective nodes can be calculated ^{59, 60}.

(6)

(7)

(8)

It is recommendable to study a "stable" network structure; a set of edges and directions that is unlikely to vary; especially for cross-sectional studies. This is because the structure of Bayesian Networks could vary when learning algorithm is performed several times. Researchers need to study stable networks whose structure is likely to be influenced by factors such as sample size ⁶¹. Bootstrapping methods can be used to evaluate the stability of network estimates in the case of partial correlation networks ⁶², and it is also possible to directly account for stability in structure learning in a similar way, where a large number of Bayesian Networks are learned from bootstrap samples and only edges that appear in a proportion of Bayesian Networks higher than some threshold are considered ^{61, 63}. In a bootstrapping method, with at least 100 bootstrap samples, edges that appear in more than 85% of networks (edge strength), and whose direction appeared in more than 50% of networks (minimum direction) are considered accurately estimated and stable enough to be included in the network and reported ^{61, 64, 65}.

Arcs, also known as edges represent the causal dependencies among variables (nodes) in the probabilistic directed acyclic graph (DAG) of a Bayesian Network structure. They connect nodes that interact and indicate the direction of influence ^{66, 67, 68}. Significant arcs can be determined using a threshold estimated from the data and referred to as a significant threshold ⁶⁹.

A systematic approach to identifying statistically significant features in a network was developed using bootstrap resampling and model averaging ^{70, 71, 72}; and summarized by ⁶⁹ as:

1. For :

a. A new data set is sampled from the data set , using bootstrap;

b. The structure of the graphical model is learned from .

2. The probability that each possible edge , is present in the true network structure is estimated as

(9)

Where is the indicator function of the event ; which implies that it is equal to 1 if , and 0, otherwise.

The observed probabilities are referred to as arc strengths, and explained as the degree of confidence that is present in the structure of the network which describes the true dependence structure of ⁷³. It is worth mentioning that the probabilities are an estimator of the expected value of the random vector which describes the presence of each possible edge in ⁷⁴. They therefore do not sum to one and are dependent on one another in a nontrivial way ⁷³, ⁷⁵.

The probability distribution of the structure of the graphical model in the space of the network structures is unknown, hence it is difficult to evaluate these empirical probabilities. The value of the confidence threshold is also is an unknown function of both the data and the structure learning algorithm. This threshold is the minimum degree of confidence for an edge to be considered significant and accepted as an edge of ⁷⁵. This leads to the use of ad hoc, pre-defined thresholds, which serve as a limitation in the identification of significant edges ⁷⁶; and a statistically motivated algorithm should be used for choosing a suitable threshold instead ⁷³.

This study adopted a solution proposed by ⁶⁹, which is implemented in bnlearn package to produce a default threshold value when a value is not specified by the analyst. They proposed that the threshold value could be computed through bootstrapping as;

If the value of the stored threshold is:

(10)

and

Then the arc strength that closely represent the unknown averaged network as

(11)

that is to say, the set of strengths that designates any arc as being significant or not with certainty. An arc is considered significant if its arc strength value is greater than the computed threshold.

Link strength (LS), is a measure of the maximum amount that evidence at the parent node of the link can affect the belief at the child, given that all the other parents of the child have some evidence ⁷⁷. It is a single number measure for each link in a Bayesian Network.

The link strength of is defined as:

(12)

where is the set of parents of and is the connection strength of .

From Figure 1, the strength of the link is the highest value of connection strength that links A to B as all the parents take on varying evidence.

(13)

A link is significant enough to be included in a DAG if the link strength is greater than 0 ^{77, 78}. This is used to determine the significance or otherwise of a link between cereal crop yield and climate change variables.

Forecasting with DBN models can be performed up to some horizon , using a sliding window ⁷⁹. Given the initial state vector with values of the variables in a system observed at instant and at many previous instances as defined by the Markovian order of the network, inference is performed to obtain the values, , that the variables are predicted to take at the next instant. All the preceding evidence is then moved forward in time. The oldest evidence is forgotten from the system and is introduced as the new evidence of the last instant to create the new state vector . The current inference step is completed at this stage, while is employed as the opening state vector of the next inference step, to predict the next state of the system. To reach the desired horizon , it is required that one should perform as many inference steps as needed. When the forecasting is completed, the values of the target variables at each instant are returned, and the mean absolute error (MAE) is calculated with the original TS if some test data are provided ^{49, 79}.

The study will employ R version 4.3.3 software from the Comprehensive R Archive Network (CRAN) project, with its corresponding r packages.

3. Results and Discussions

The DBN network structure as learned is shown below in Figure 2. It shows the network connection between climate variables and crop yield at the previous time slides. It also indicates a connection from the previous time slide to the current time slides, showing how previous and current time slides influence cereal crop production.

The structure indicates by the direction of the arrows that there is an association within the time slides, where climate variables in the previous slides show causal relationships with cereal crop yield at that time slide (Cropc_t_1). The side with the blue balls represents the current time slide conditioned on previous time slides. There is a connection that links all the climate variables at previous time slides to the crop yield at the current time slide, including cereal crop yield in previous time slides. There is also an indication of climate variables at the current time slide (intra-slide) having links with crop yield. These links are all represented by directed arcs toward cereal crop yield, which indicate that the variables at the tail have some effect on those at the head. Climate variables at each time slide have an effect on cereal crop yield at the current time slides from our data, according to the model as fit. That is, cereal crop yield is associated with the previous and current climate variables considered.

Upon learning the structure of the DBN, the parameter learning is done and the estimates are presented in Table 1 below. The table indicates the variables at the previous and current time slides with their corresponding partial correlations and link strengths.

From the Table 1, we observe the values of the various nodes at the different time slides under the partial correlation column which indicates how these nodes can be used as predictors of crop yield. The odds ratio column represents the classical odds ratio (OR), and it is a measure of the relationship between the random variables at the different time slides. It indicates that all the variables from the previous time slide are positively associated with the current crop yield with values greater than one, except windspeed which has a value less than one, and hence negatively associated with the current crop yield. The OR values for the intra-slide relationship at the current time slide indicate that all the climate values have a positive relationship with crop yield at the current time-slide, except temperature (temp_t_0) which has a negative relationship.

The link strength column shows values that represent the measure of the highest amount that evidence at the parent node of a link can effect the belief, given that every one of the other parents of the child has some evidence. It indicates the significance of the link if this value is either greater than or more than 0.5. From Table 1 therefore, the link strength indicates significant connections or influence of the climate variables at the different time slides with crop yield (Cropc_t_0) at the current time slide except for crop yield (Cropc_t_1) and windspeed (Windsc_t_1) at the previous time slides.

With the intra-slide interaction at the previous time slide, we observe that all the climate variables are positively and significantly associated with the cereal crop yield, except humidity (Humc_t_0) which is positively but not significantly associated with crop yield.

Figure 3 represents the DAG indicating the stability of the DBN model fitted. With a ten thousand resampling process, the DAG indicates that the arcs and arc directions remain stable. It indicates that in the face of increased data points and additional variables, the model will still fit well. Thick arcs in the diagram indicate that there is a negative and non-significant relationship.

The purpose is to forecast the number of cereals in metric tons to be produced for a specified year, since obtaining accurate predictions will help enhance production efficiency decision-making ahead of time before and during cultivation. The model indicates a mean absolute error of 0.0429 which stabilizes towards zero as the number of predicted instances increases. Figure 4 shows predicted values (red) almost converging with true values with time.

4. Conclusions

The objective of the study is to investigate and model the effect of climate change on cereal crop yield in Northern Ghana, using the Dynamic Bayesian Network model. Using the research data (cereal crop yield and climate variables) for a thirty-one-year period, the Dynamic Bayesian Networks model was fit and model diagnostics were assessed.

The results indicated the existence of intra and inter-slide connections in the model. It was determined that climate variables at each time slide had shown to have some relationship with cereal crop yield at the current time slide. That is, cereal crop yield is associated with the previous and current climate variables. It also revealed that all the climate variables were positively and significantly associated with the cereal crop yield, except humidity (Humc_t_0) which was positively but not significantly associated with crop yield. Model diagnostics processes indicated that the model fits well with the data. The model also has a Mean Absolute Error (MAE) value of 0.0429 which indicates better forecasting accuracy and has the tendency to perform well in understanding or assessing the relationship between cereal crop yield and climate change variables in Northern Ghana.

The model results showed that climate variables are significantly linked with yield at the intra and inter-slide levels. It was also proven to be sensitive enough to capture association even with small data and showed stability to out-of-sample variations. It was further revealed that there exist significantly positive relationships between cereal crop yield and climate change variables. The results imply that, with data on current climate variables, we can predict the next cereal crop yield with reduced uncertainty and better forecasting ability while holding other factors and inputs constant. Again, it was established that any misspecification in rainfall, temperature, and humidity has a greater potential effect on yield since cereal crop yield has shown to be sensitive to these.

The study therefore recommends the use of this approach in understanding the relationship between cereal crop yield and climate change. Also, climate change impacts should always be taken into consideration and factored into decision-making processes before and during cultivation; since this study has shown that such relationships exist. Besides, since other factors could influence cereal crop yield, it is recommended that the study be extended to include other variables such as vegetative cover, area cultivated, chemical use, etc. when data becomes available.

References