# Linking Theories and Methods in Cognitive Sciences via Joint Embedding of
the Scientific Literature: The Example of Cognitive Control

Morteza Ansarinia
University of Luxembourg
Department of Behavioural and Cognitive Sciences
<EMAIL_ADDRESS>Paul Schrater
University of Minnesota, Twin Cities
Depts. of Computer Science and Psychology
<EMAIL_ADDRESS>Pedro Cardoso-Leite
University of Luxembourg
Department of Behavioural and Cognitive Sciences
<EMAIL_ADDRESS>
###### Abstract

Traditionally, theory and practice of Cognitive Control are linked via
literature reviews by human domain experts. This approach, however, is
inadequate to track the ever-growing literature. It may also be biased, and
yield redundancies and confusion.

Here we present an alternative approach. We performed automated text analyses
on a large body of scientific texts to create a joint representation of tasks
and constructs. More specifically, 385,705 scientific abstracts were first
mapped into an embedding space using a transformers-based language model.
Document embeddings were then used to identify a task-construct graph
embedding that grounds constructs on tasks and supports nuanced meaning of the
constructs by taking advantage of constrained random walks in the graph.

This joint task-construct graph embedding, can be queried to generate task
batteries targeting specific constructs, may reveal knowledge gaps in the
literature, and inspire new tasks and novel hypotheses.

Keywords: Cognitive Control; Natural Language Processing; Cognitive
Constructs; Cognitive Tasks;

## Introduction

A key challenge in cognitive sciences, and in particular cognitive psychology
and neuroscience, is to make sense of observable phenomena (i.e., behavior) in
terms of theoretical constructs. Consider for instance Cognitive Control
(CC)—a broad construct that comprises many components and engages multiple
mechanisms which collectively aim to describe goal-directed behavior in a
complex, uncertain world.

CC is a major construct in cognitive sciences: In the year 2021 alone, PubMed
indexed 974 papers with the term ”Cognitive Control” in the title or
abstract—an average of 3 papers per day. To understand CC, researchers have
introduced a variety of theoretical constructs and conceived numerous
cognitive tasks (see @baggetta2016). However, the relationships between and
within related constructs and tasks are not always clear. For example, because
they are ”measured” using the same set of tasks (e.g., Stroop, N-back, Digit
Span, Stop-Signal, Task Switching), it seems reasonable to assume that
cognitive control Botvinick Cohen (2014), executive functions Baggetta
Alexander (2016), attentional control Rey-Mermet . (2021), and self-regulation
Enkavi . (2019) are somewhat equivalent constructs; yet, they are not widely
considered equal Nigg (2016).

Traditionally, the meaning and relationships between constructs and tasks are
conceptualized in extensive literature reviews conducted by human experts. In
this approach, researchers ”manually” read, synthesize, and criticize the
literature and write reviews or reports describing their understanding.
Following such reviews, CC is viewed as interactions between generic core
processes (e.g., inhibition, flexibility, working memory, and interference
control in @diamond2013), interactive componential Badre (2011), tasks-
specific processes driven by goals Logan (2017); Doebel (2020), or optimal
parameterization of naturalistic tasks Botvinick Cohen (2014). This approach
has been invaluable but it may also yield biased results Brick . (2021); Beam
. (2021) and seems inadequate to track the ever-growing literature and stay
current. In this context, modern machine learning methods may provide useful
and complementary insights.

When considering terms in the literature, there are two major impediments to
creating consistent construct-task associations: construct hypernomy when
conceptualizing CC and task impurity when operationalizing it. Construct
hypernomy occurs when description of the same construct varies across
different contexts due to the way it is assessed. It creates different
meanings of the same concept. ”Attentional Control”, for example, likely means
something different in @ahissar1993 (as measured by low-level perceptual
tasks) than it does in @burgoyne2020 (as measured by complex cognitive tasks).
Task impurity, on the other hand, refers to the idea that performance on a
task loads onto multiple constructs (i.e., there is not a one-to-one mapping
between constructs and tasks). Because of the impurity, no task taps into just
one isolated construct. Performance in the Backward Digit Span, for instance,
involves short-term memory, visual perception, sustained attention and working
memory, to name just a few. The consequence is that constructs lack a
consistent, groundable semantic content, corrupting interpretations of neural
and cognitive research that depend on them.

Construct hypernomy and task impurity are quite common in CC research because
complex concepts like cognitive control manifest themselves differently across
different individuals and contexts Burgoyne Engle (2020). For that,
researchers often use multiple tasks in their studies and apply statistical
methods such as latent factor analysis to discern underlying constructs.
Nevertheless, the resulting latent models of CC are rarely agreed upon, as is
the selection of tasks (see, for example, @reymermet2021; @doebel2020,
@enkavi2019, @nigg2016).

Ambiguous associations of constructs and tasks make it hard to interpret past
results, hinder scientific progress and the development of effective
interventions. With the advent of scalable machine learning, however,
construct-task associations may be clarified. The goal of this paper is to
approach the conceptual richness of a large body of scientific works and take
advantage of recent context-aware language models in machine learning to
clarify the association of CC tasks and constructs. More specifically, we
collect and analyze scientific texts about CC tasks and constructs and encode
text data into rich semantic embeddings using transfer learning. Transfer
learning exploits the rich representations generated by natural language
models trained to faithfully represent contextual meaning—unlike traditional
bag-of-word or clustering techniques. Similarities between embedded
representations are then used to build up a hypergraph Battiston . (2021) that
connects tasks and constructs.

First, we show that this hypergraph representation regrounds constructs on
tasks and provides nuanced meaning of the constructs, ultimately demonstrating
construct hypernomy. Second we show that pulling theoretical and experimental
literature into overlapping components of a hypergraph may greatly benefit
researchers: the joint task-construct embeddings can be queried to generate
special-purpose task batteries, it may reveal knowledge gaps, inspire the
design of new experiments and yield novel hypotheses regarding the structure
and function of CC. This empirical and descriptive model of the literature,
rather than expert-driven ones, may also be used in future applications to
enhance knowledge searches (see @beam2021 for a comparison of a data-driven
mapping of the literature and expert-driven knowledge frameworks like DSM and
RDoC).

## Methods

Data. We created a lexicon of CC-related terms (172 terms, of which 72 were
task names and 100 were construct names) based on the previously published
work on Cognitive Control Barch . (2009), Attentional Control Bastian . (),
Executive Functions Baggetta Alexander (2016); Diamond (2013), and Self
Regulation Enkavi . (2019). Each term in the lexion was associated with a
PubMed-specific search query by which papers with the term in their title or
abstract were retrieved. This resulted in a dataset of loosely labeled
documents, each labeled by one or more lexicon terms (n=522,972 hits, of which
385,705 were unique). For the purpose of the current analyses we only retained
the title and abstract of the papers, along with the lexicon terms that were
used to retrieve them. Having multiple labels per document was crucial to
quantify the co-appearance of the terms in the literature. After the documents
were collected we removed 14 terms from the lexicon because they yielded too
few documents to support cross validation splits ($n<5$).

Analysis. To understand the relationships among and between tasks and
constructs, our goal is to build graphs that represent tasks and constructs as
nodes and measure similarity/distance between them as edges. Graph $G$ can be
used to jointly infer embeddings of both construct and task nodes in a shared
vector space, in that relative closeness of two nodes is estimated by the
similarities of node attributes as well as the shared neighbors in the graph.
Heterogeneous graph $G=(V_{tasks}\bigcup V_{constructs},E)$ is defined by its
two types of nodes, $V_{tasks}$ and $V_{constructs}$, labeled by either a task
or a construct term, while the weighted edges, $E$, represent the links
between two or more nodes, reflecting similarity of the corresponding terms in
the literature. Node attributes being relevant scientific texts, the existence
and weight of a link between two nodes is predicted by the similarity of
corresponding node attributes; the higher the similarity between node
attributes, the higher the chance of the nodes being associated. The core
problem becomes learning task and construct attribute embeddings that predict
co-occurrence and semantic similarity measures. We used the following steps to
create the graph $G$ from the collected scientific texts.

(A)

(a) *

(B)

(b) *

(C)

(c) *

(D)

(d) *

Figure 1: (A) Introducing new tasks (task innovation) and constructs (concept
innovation) is characterized by a burst followed by declining innovation. (B)
Task and Construct occurrences in publication abstracts are temporally
decoupled. Time to operationalize constructs (blue) is the time between the
first occurrence of a construct and the first co-occurrence of that construct
with any tasks, while Time to conceptualize tasks (orange) is the time between
the first occurrence of a task and the first co-occurrence of that task with
any of the construct. (C) The majority of the literature only used one task in
their studies, showing a lack of multitask design of experiments. (D) While
the number of papers published each year increases exponentially, the number
of tasks per study remains fairly constant across time. Figure 2: Task-
Construct hypergraph: representations of control-related constructs as
hyperedges (vertical black lines) over a subset of tasks (nodes). Construct
hypernomy is reflected as overlapping hyperedges (e.g., green regions), and
task impurity as nodes scattered over multiple hyperedges (e.g., blue region).
Distances between nodes are not meaningful. Nodes are reorganized for visual
clarity and only a subset of the graph is displayed.

The data collection resulted in a dataset of 385,705 unique, but loosely-
labeled, abstract texts, all of which were then encoded into embeddings of
1024 dimensions using a pre-trained transformer language model (GPT-3 Ada for
text similarity embedding; see @brown2020). The language model transformed raw
texts into 1024-dimensional vectors, gpt3-embedding, representing semantic
similarity between two or more pieces of text. Since keeping the original
structure of the text was important for the model to understand the context,
we did not preprocess the raw text. To convert text similarity into a shared
topic representation (which improves relating task and construct text
embeddings), we applied Top2Vec topic modeling Angelov (2020) to the
gpt3-embedding which projected them into a space of 473 dimensions, i.e.,
topic-embeddings. Each column of the topic-embedding matrix represents a
topic, and element $ij$ shows the probability of assigning document $i$ to the
topic $j$. Realigning the gpt3-embedding into topic-embeddings improved the
quality of the dataset for a number of reasons. First, it improves the quality
of the labels in the dataset by discarding outlier documents. These are
documents that belong to no topics of interest or are assigned to irrelevant
topics (e.g., genetics)—after removing outlier documents, 293’014 unique
documents remained for further analysis. Second, topic modeling allows one to
extract a useful, interpretable representation of the documents, as each
dimension of the topic-embedding shows the probability of assigning a document
to a topic while being faithful to the contextual representation of the
documents in the gpt3-embedding space. This generates a digraph between nodes
representing lexicon terms and the topic-embedding vectors.

To convert this into a construct-task graph, we grouped lexical terms
associated with construct and tasks to generate graph nodes. To compute topic-
similarity between groups of lexical terms associated with each construct or
task node, we fitted a multivariate normal distribution over the topic vectors
of each node separately and then calculated the distance between all nodes as
measured by the Jensen-Shannon divergence of those node-level distributions.
This step added edges to the graph, G, with edges weighted by the inverse
distance of nodes in the JS-divergence matrix.

To learn a representation of the graph that only preserves paths from tasks to
constructs and vice versa, we then applied Metapath2Vec (1000 random walks of
step size 100, accompanied by skip-gram Word2Vec embedding of size 128 and
maximum window size of 5; as recommended in @ruch2020). The Metapath2Vec
embedding encodes random walks of specific patterns in a heterogeneous graph,
here patterns being alternating random walks between task and construct nodes.

Finally, by applying HDBSCAN soft clustering to the node attributes and
thresholding the edges (discarding all the edges weighed within one standard
deviation from median), we transform the graph G to a homogeneous hypergraph,
i.e., nodes are now only of type task, while constructs are hyperedges that
group a subset of tasks in overlapping clusters.

## Results

We used a variety of data-driven approaches to collect and understand CC
publications. Briefly, we a) created an all-inclusive lexion of construct and
task terms, b) queried PubMed to collect relevant abstract texts, c)
vectorized all the raw texts using GPT-3 Sentence Similarity Embedding, an
unsupervised pre-trained language model, d) applied Top2Vec topic modeling
technique to all the document embeddings together and identified dimensions of
a useful latent space, i.e., topics. We then created a graphical
representation of the lexicon terms, i.e., task-construct graph, and used them
to predict the association between terms.

The richness of tasks and constructs in the literature. Although there are
many task and construct terms, their relative frequencies differ widely. For
example, ”Stroop Task” is mentioned 8,003 times in the period 1973–2022 while
”Delay Discounting Task” was only mentioned 466 times over the same period of
time. The use of each term tends to increase over time. Interestingly the rate
at which new constructs and tasks are introduced does not follow the same
curve as the number of publications in the field; rather there seems to have
been a peak of innovation for constructs around 1980 and for tasks around the
year 2000 (Figure 1A). Such patterns, visible in simple descriptive statistics
(Figure 1), may provide interesting insights into understanding the maturity
and vitality of a research field.

Regrounding constructs on tasks. It took on average 7 years for the constructs
to be explicitly associated with a task (see Figure 1B). The meaning of a
theoretical construct may change across time and gain clarity and precision
with new empirical measures and cognitive tasks being used by the research
community to flesh out the construct. A core idea in this paper is that by
evaluating how constructs are operationalized (i.e., linked to cognitive
tasks) key insight can be gained about what a construct means. Grounding the
definition of constructs on tasks provides a nuanced meaning of constructs
that relies on observable measures. It also allows the computation of useful
measures on constructs (e.g., specificity) and on between pairs of constructs
(e.g., measures of redundancy, similarity, and distance). To investigate the
relationships between cognitive constructs, we use hyperedges in the task-
construct graph as a measure of similarity, indicating the extent to which a
construct hyperedge can be reconstructed by neighboring tasks.

Construct hypernomy. The task-construct graph readily demonstrates construct
hypernomy and task impurity in the CC literature. We first sought hypernomy as
highly overlapping hyperedges of seemingly incompatible constructs, as well as
a high degree of task nodes with neighboring constructs as a measure of the
task impurity. Figure 2 illustrates overlapping hyperedges of the most popular
constructs where hyperedges for Cognitive Control, Executive Control,
Behavioral Control, Central Executive, and Attentional Control are overlapping
and identical.

Figure 3: Associations between tasks and constructs minimally overlap across
scientific disciplines. Rose plots show the relative association between
constructs and tasks, with each color representing a different field. Lack of
overlap between the ”spikes” indicates disjoint operationalizations across
fields. Figure 4: Pairwise distances between the 25 most popular cognitive
control tasks as measured by the symmetric Jensen-Shannon divergence of two
multivariate normal distributions of their node attributes in the task-
construct graph. Higher divergence indicates higher dissimilarity between
corresponding scientific texts. Task-task distances may for example provide a
data-driven proxy for predicting and explaining transfer effects in cognitive
training research.

Task inconsistency across disciplines. A major source of hypernomy stems from
descriptions and measurements of the constructs often being inconsistent
across scientific communities. To test this idea, we sought to determine
whether construct hyperedges, and their task associations, vary across four
cognitive disciplines (psychology, neuroscience, cognitive science, and social
science). Using the same method described in the analysis section we created
four discipline-specific graph embeddings. The only difference was that
publications were grouped by discipline, which was determined by searching for
the terms ”social”, ”psycho”, ”neur”, or ”cognit” in the journal titles.
Constructs that have inconsistent task associations across the disciplines are
hypernomic (Figure 3).

Refactoring tasks and constructs. Designing effective assessments of CC can be
challenging for a number of reasons. Participants have limited time to spend
on cognitive tasks. 1) If these tasks are poorly selected, performance on
these tasks may not be very informative (e.g., measures are conceptually
redundant); 2) If only one task is used, the inferential resolution of
performance to construct is very limited. Thus in order to be able to make
specific theoretical claims about CC it is necessary to use multiple, well-
chosen tasks in experiments. This is currently not the case. As shown in
Figure 1C, most research uses only one task. In fact, only 17 percent of
publications used 2 or more tasks. The task-construct graph presented here may
facilitate novel experimental designs of such multi-task, max-information
experiments by providing a similarity-based space in which tasks can be
identified, and grouped, by the overlapping subgraphs (i.e., constructs) that
they belong to.

In the task-construct graph, two tasks are similar if they share identical
neighbors, i.e., constructs. And tasks cover a set of constructs if their
union set overlaps the corresponding hyperedges of the constructs. These
principles equip researchers with sound and quantified methods to refactor
tasks (e.g., discard redundant tasks, quantitatively measuring similarity of
tasks via constructs, and performing set operations on a group of tasks). Such
a refactored set of tasks controls the construct-redundancy of tasks and will
shorten the time required to complete comprehensive assessments. It provides a
method to design a task battery to effectively cover constructs (i.e., minimal
redundancy while measuring different facets of the constructs).

Sparsity in the task space. There are numerous cognitive tasks in the
literature; how these tasks relate to each other remains unclear. There are
many cognitive control tasks that are rarely used (see @baggetta2016), and
even fewer used in combination with other tasks. Even when tasks were used
together, their relationship might still be unclear. The question of how tasks
relate to each other is key in the cognitive training domain where researchers
aim to train cognitive abilities in general rather than performance on a
specific task. In that context, a common point of disagreement is to predict
and interpret transfer effects (i.e., how much training in task A improves
performance in task B). A measure of distance between tasks based on their
grounding on constructs may provide an objective foundation to understand
these transfer effects—the task-construct graph embedding proposed here
provides a means to compute such inter-task distances.

To quantify the distance between two cognitive tasks, we compute the Jensen-
Shannon divergence between their node embeddings in the task-construct graph.
Figure 4 shows, for example, that the Trail Making Task is relatively close to
the Digit Span Task, suggesting its training effects transfer more easily to
the Digit Span Task than to tasks such as the Discounting Task.

Distance between the task nodes can also allow us to identify gaps in the task
space: gaps may be visible as disconnected graph components. Identifying such
gaps may reveal opportunities to develop new useful tasks. Alternatively,
there may only exist associations between groups of tasks and groups of
constructs—i.e. the task-construct associations are not atomic. This reflects
a lack of purity in the tasks or constructs or both that might be improved by
refactoring constructs or decomposing tasks into components.

Querying the graph embedding for task batteries targeting specific cognitive
constructs. Some studies use batteries of tasks that together address a
research question and measure one or more constructs from several viewpoints.
The process of building such task batteries can be facilitated by leveraging
the task-construct graph embedding; one can query the graph for an array of
tasks spanning a given set of constructs. The joint embedding translates
queries into arithmetic operations in the embedding space (positive samples
and negative samples), allowing for more explicit and visible decisions.

Query operations on the task-construct graph are made possible by using the
underlying node embedding vectors extracted as a part of Metapath2Vec graph
embedding. Queries include, for example, prioritizing tasks for a given
construct, or a set of tasks for a set of constructs. To prioritize tasks for
a construct, the task-construct graph looks for task nodes that are closest to
the simple mean of the queried construct, e.g., in terms of sum of weighted
node embeddings. And for a list of tasks for multiple constructs, find the
minimum spanning tree that covers all the queried construct hyperedges. For
example, if one queries (Reward Processing + ReversalLearning - GoNoGo -
SortingTask), one will get the recommendation to use the BART, GiftDelay,
BalanceBeam Baggetta Alexander (2016), and StimSSS Enkavi . (2019) tasks,
which are ordered by the cosine similarity between the mean vector of the
query and the task vectors in the graph embedding model.

## Implications

Ambiguous meanings and relationships between cognitive tasks and constructs
call for a more rigorous way to handle constructs—an obvious solution would be
to adopt a more formal notation and refer to specific knowledge models (e.g.,
ontologies). The knowledge model must be flexible enough to capture a wide
range of association between constructs and tasks. The proposed task-construct
graph embedding provides a useful representation of the cognitive control
literature built upon topic embedding. In this representation, association of
two entities, e.g., task-construct, relies on shared topics as well as the
walks between them in a graph representation. By predicting links using topic
embeddings of the nodes, we find most similar aspects of, for example, two
constructs, a similarity that could be explainable in natural language.

A consistent, sound, and parsimonious framework of CC has been desired from
the beginning. Yet, the growing number of publications and newly introduced
constructs makes it impossible to integrate them into a bigger picture. While
researchers may disagree on theoretical perspectives and thus on which terms
to use, they generally might agree on the fact that if two constructs are
”measured” by the same tasks, the constructs must be somewhat related. We
proposed a joint embedding of constructs and tasks (based on scientific texts
in a graph representation) to drive a more nuanced interpretation of the
constructs by regrounding abstract constructs on the concrete set of
observable tasks.

The proposed graph-based embedding enables explanatory reasoning driven by
scientific texts. Unlike expert-driven models, the models reason regardless of
the preferences in research; yet it is not clear whether other kinds of biases
are addressed as the knowledge source and pre-trained language model are
themselves produced by humans. By scaling up the knowledge model to a large
body of available texts, the model is able to encapsulate even more aspects of
Cognitive Control, and in general, multidisciplinary research.

Disagreements about the meaning of a construct are partly explained by
differences in how we interpret responses to a particular task. By focusing on
the co-occurrence of task and construct names in scientific texts, our
approach implicitly makes strong assumptions about the relationship between
abstract constructs and their imperfect but observable measures. The
limitations of the present work can be partially addressed by expanding the
hypergraph to include, for example, concepts such as brain mechanisms,
research communities, and analysis techniques.

Explainable symbolic AI and machine learning have been long in debate to model
knowledge. Regardless of the specific topic discussed here (i.e., Cognitive
Control), the proposed model can be seen as an effort to connect symbolic
modeling (as in ontologies) and machine learning (as in embeddings). Our
method informs an ontology of scientific texts using context-aware embeddings
that are extracted from a loosely-labeled body of scientific texts requiring
minimal human input. It is an automated pipeline that only requires a lexicon,
builds on large-scale language models and that can scale to millions of
documents, making it a viable approach to meaningfully monitor the scientific
literature continuously and extensively.

## Acknowledgments

This research was supported by the Luxembourg National Research Fund
(ATTRACT/2016/ID/11242114/DIGILEARN) and (INTER
Mobility/2017-2/ID/11765868/ULALA).

## References

  * Ahissar Hochstein (1993) @ahissar1993Ahissar, M. Hochstein, S.  1993\.  Attentional control of early perceptual learning Attentional control of early perceptual learning.  Proceedings of the National Academy of Sciences9012.
  * Angelov (2020) @angelov2020Angelov, D.  2020\.  Top2Vec: Distributed Representations of Topics Top2Vec: Distributed Representations of Topics.  arXiv.
  * Badre (2011) @badre2011Badre, D.  2011\.  Defining an Ontology of Cognitive Control Requires Attention to Component Interactions Defining an Ontology of Cognitive Control Requires Attention to Component Interactions.  Topics in Cognitive Science32.
  * Baggetta Alexander (2016) @baggetta2016Baggetta, P. Alexander, PA.  2016\.  Conceptualization and Operationalization of Executive Function Conceptualization and Operationalization of Executive Function.  Mind, Brain, and Education101.
  * Barch . (2009) @barch2009Barch, DM., Braver, TS., Carter, CS., Poldrack, RA. Robbins, TW.  2009\.  CNTRICS Final Task Selection: Executive Control CNTRICS Final Task Selection: Executive Control.  Schizophrenia Bulletin351.
  * Bastian . () @bastian2021Bastian, CCv., Blais, C., Brewer, GA., Gyurkovics, M., Hedge, C., Kałamała, P.Wiemers, EA.  .  Advancing the understanding of individual differences in attentional control: Theoretical, methodological, and analytical considerations Advancing the understanding of individual differences in attentional control: Theoretical, methodological, and analytical considerations.  PsyArXiv.  10.31234/osf.io/x3b9k
  * Battiston . (2021) @battiston2021Battiston, F., Amico, E., Barrat, A., Bianconi, G., Arruda, GFd., Franceschiello, B.Petri, G.  2021\.  The physics of higher-order interactions in complex systems The physics of higher-order interactions in complex systems.  Nature Physics1710.
  * Beam . (2021) @beam2021Beam, E., Potts, C., Poldrack, RA. Etkin, A.  2021\.  A data-driven framework for mapping domains of human neurobiology A data-driven framework for mapping domains of human neurobiology.  Nature Neuroscience.
  * Botvinick Cohen (2014) @botvinick2014Botvinick, MM. Cohen, JD.  2014\.  The Computational and Neural Basis of Cognitive Control: Charted Territory and New Frontiers The Computational and Neural Basis of Cognitive Control: Charted Territory and New Frontiers.  Cognitive Science386.
  * Brick . (2021) @brick2021Brick, C., Hood, B., Ekroll, V. de Wit, L.  2021\.  Illusory Essences: A Bias Holding Back Theorizing in Psychological Science Illusory Essences: A Bias Holding Back Theorizing in Psychological Science.  Perspectives on Psychological Science.
  * Brown . (2020) @brown2020Brown, TB., Mann, B., Ryder, N., Subbiah, M., Kaplan, J., Dhariwal, P.Amodei, D.  2020\.  Language Models are Few-Shot Learners Language Models are Few-Shot Learners.  arXiv.
  * Burgoyne Engle (2020) @burgoyne2020Burgoyne, AP. Engle, RW.  2020\.  Attention Control: A Cornerstone of Higher-Order Cognition Attention Control: A Cornerstone of Higher-Order Cognition.  Current Directions in Psychological Science296.
  * Diamond (2013) @diamond2013Diamond, A.  2013\.  Executive Functions Executive Functions.  Annual Review of Psychology641.
  * Doebel (2020) @doebel2020Doebel, S.  2020\.  Rethinking Executive Function and Its Development Rethinking Executive Function and Its Development.  Perspectives on Psychological Science154.
  * Enkavi . (2019) @enkavi2019Enkavi, AZ., Eisenberg, IW., Bissett, PG., Mazza, GL., MacKinnon, DP., Marsch, LA. Poldrack, RA.  2019\.  Large-scale analysis of test–retest reliabilities of self-regulation measures Large-scale analysis of test–retest reliabilities of self-regulation measures.  Proceedings of the National Academy of Sciences11612.
  * Logan (2017) @logan2017Logan, GD.  2017\.  Taking Control of Cognition: An Instance Perspective on Acts of Control Taking Control of Cognition: An Instance Perspective on Acts of Control.  American Psychologist729.
  * Nigg (2016) @nigg2016Nigg, JT.  2016\.  Annual Research Review: On the relations among self-regulation, self-control, executive functioning, effortful control, cognitive control, impulsivity, risk-taking, and inhibition for developmental psychopathology Annual Research Review: On the relations among self-regulation, self-control, executive functioning, effortful control, cognitive control, impulsivity, risk-taking, and inhibition for developmental psychopathology.  Journal of Child Psychology and Psychiatry584.
  * Rey-Mermet . (2021) @reymermet2021Rey-Mermet, A., Singmann, H. Oberauer, K.  20213\.  Neither measurement error nor speed-accuracy trade-offs explain the difficulty of establishing attentional control as a psychometric construct: Evidence from a latent-variable analysis using diffusion modeling Neither measurement error nor speed-accuracy trade-offs explain the difficulty of establishing attentional control as a psychometric construct: Evidence from a latent-variable analysis using diffusion modeling.  PsyArXiv.  https://doi.org/10.31234/osf.io/3h26y
  * Ruch (2020) @ruch2020Ruch, A.  2020\.  Can x2vec save lives? Integrating graph and language embeddings for automatic mental health classification Can x2vec save lives? Integrating graph and language embeddings for automatic mental health classification.  Journal of Physics: Complexity13.


# Crop yield probability density forecasting via quantile random forest and
Epanechnikov Kernel function

Samuel Asante Gyamerah  Philip Ngare  Dennis Ikpe  Pan African University,
Institute for Basic Sciences, Technology, and Innovation, Kenya  University of
Nairobi, Kenya  Michigan State University, USA

###### Abstract

A reliable and accurate forecasting model for crop yields is of crucial
importance for efficient decision-making process in the agricultural sector.
However, due to weather extremes and uncertainties, most forecasting models
for crop yield are not reliable and accurate. For measuring the uncertainty
and obtaining further information of future crop yields, a probability density
forecasting model based on quantile random forest and Epanechnikov kernel
function (QRF-SJ) is proposed. The nonlinear structure of random forest is
applied to change the quantile regression model for building the probabilistic
forecasting model. Epanechnikov kernel function and solve-the equation plug-in
approach of Sheather and Jones are used in the kernel density estimation. A
case study using the annual crop yield of groundnut and millet in Ghana is
presented to illustrate the efficiency and robustness of the proposed
technique. The values of the prediction interval coverage probability and
prediction interval normalized average width for the two crops show that the
constructed prediction intervals capture the observed yields with high
coverage probability. The probability density curves show that QRF-SJ method
has a very high ability to forecast quality prediction intervals with a higher
coverage probability. The feature importance gave a score of the importance of
each weather variable in building the quantile regression forest model. The
farmer and other stakeholders are able to realize the specific weather
variable that affect the yield of a selected crop through feature importance.
The proposed method and its application on crop yield dataset are the first of
its kind in literature.

###### keywords:

climate change , crop yield uncertainty , crop yield forecasting , quantile
random forest , kernel density estimation , Epanechnikov kernel

label1label1footnotetext: correspondence<EMAIL_ADDRESS>
##  1 Introduction

The agriculture sector is seen as one of the biggest emitter of greenhouse
gases and concurrently a major sector that is affected by climate change. In
reviewing the factors that affect crop growth, productivity, and yield, [1]
indicated that soil moisture, the availability of soil nutrients, and solar
radiation are the top three factors that limit the growth of crops and hence
limit the yield of crops. Changes in the surface temperature, humidity, and
rainfall affects the moisture content of the soil and the level of nutrients
in the soil [2]. Hence, there is a direct effect of climate change on crop
growth, productivity and yield [3, 4, 5]. This has significantly affected the
yield of most crops causing economic and food security risks in most
developing and under-developed countries [6, 7].
Small-scale farmers in most developing and under-developed countries cannot
mitigate the effects of climate variability [6]. This mostly have great impact
on the farming activities and economy of these farmers. For this reason, an
effective and reliable insurance is needed to hedge farmers and stakeholders
from the peril of weather uncertainties. Traditional insurance for
agricultural risk management is not patronized in most developing countries
because of high premiums, loss adjustments, moral hazards, adverse selections,
and complex information requirements [8]. However, weather derivatives and
index-based insurance such as area-yield and weather index insurance are seen
as effective risk management tools in the agricultural sector for both small-
and large-scale farmers in developing/under-developed countries. Accurate
forecasting of crop yields is a principal component for the ratemaking process
in the derivative and index-insurance markets. An accurate and mathematically
tractable crop yield forecasting model is important for the farmer,
policymakers, the government, field managers, and industry players in decision
making process [9].
To improve the performance of the methods used in crop yield forecasting, a
lot of research have been done in recent decades. A number of literatures
based on statistical models have been used to predict the yield of crops [10,
11, 12]. These literatures have served as a substitute to process-based
models, which always involve a comprehensive data on the conditions of the
soil, cultivar, and management. [13] used different time series models (simple
and double exponential smoothing, Damped-Trend Linear Exponential Smoothing
and autoregressive moving averages (ARMA)) to predict maize yield in five
communities in Ghana. The authors concluded that the ARMA model was more
robust than the other time-series models. However, time series models such as
moving averages, simple and double exponential smoothing, quadratic and linear
regression perform poorly in predicting crop yields [14, 15]. These
statistical predictions suffer from different sources of error like variations
in weather variables. In most statistical methods, there are always little or
no interaction between the weather variables used for prediction. However,
crop yield and weather variables are highly nonlinear and there are
interactions between weather variables [16, 17]. Hence, using statistical
methods will be computationally costly and may not lead to an optimal
performance especially when there are cases of extreme events. Alternative to
statistical models, are the emergence of machine learning (ML) techniques such
as random forest, support vector regression, and neural networks. These ML
techniques are able to capture the nonlinearity of crop yields and its
predictors [18, 19, 20]. Hybrid methods between statistical and machine
learning approach are seen to improve the accuracy of predictions in most
forecasting problems [21, 22].
Generally, most statistical and machine learning literatures on crop yield
forecasting are based on point forecasts [11, 13, 18, 19]. Point forecasts
give an estimate of the future crop yield for each time horizon and do not
convey any information about the uncertainty of the predictions. Alternative
to point forecasting is prediction interval. Prediction interval builds a well
calibrated upper and lower bound between which the future unknown prediction
value is expected to lie with a certain prescribed probability [23]. Due to
the increased uncertainty of weather in recent years as a result of climate
change [24, 25], point and interval predictions are not able to predict the
yield of crops accurately. Unlike point prediction and prediction interval,
probability density forecasting can measure the uncertainty of crop yields by
building probability density function of the forecasting results. With the
growing increase in climate change and its effects on crop yields, it is
imperative to measure the uncertainties associated with crop yield forecast.
In probability density forecasting, Kernel density estimation (KDE) is very
important in the density estimation process. KDE is a non-parametric method of
estimating the distribution of a dataset without prior assumptions of the
datasets [26]. Appropriate choice of bandwidth for a kernel density estimator
is of crucial importance to the density function of random variables. Quantile
regression (QR) can be used to construct a nonparametric probability density
forecast. Given one or more covariates, QR generalizes the theory of a
univariate quantile to a conditional quantile. Because of the robustness of QR
in handling outliers in explained measurements, it is widely used for
regression analysis in the areas of econometrics and statistics [27].
Conventional linear QR is however unable to deal with complex non-linear
problems [28]. To explore non-linear functions for QR, [29] proposed a
quantile random forest (QRF) model, which combines the advantages of random
forest and quantile regression models. In furtherance to the application of
Meinshausen proposed model, [22] proposed a hybrid semi-parametric quantile
regression forest to estimate the non-linear relationship in multi-period
value-at-risk. They concluded that their proposed method was more accurate
compared to common distributions like normal distribution. In the area of
medicine, [21] applied quantile regression forest to Cancer Cell Line
Encyclopedia (CCLE) dataset to give a point and interval prediction. The
quantile regression forest improved the accuracy of prediction of drug
response. However, these literatures [e.g. 21] are not able to obtain the
probability density functions of the future response variable in a single QRF
model. Hence the need to combine quantile regression forest and kernel density
estimation to estimate the probability density functions of future crop
yields. To obtain a complete crop yield probability density curve,
Epanechnikov kernel function and solve-the-equation plug-in approach of
Sheather and Jones (SJ) bandwidth selection method [30] are combined with
quantile random forest (QRF) model. Our proposed method (QRF-SJ) will help in
obtaining a complete conditional probability density for different time
horizons by selecting a suitable bandwidth and kernel function.
To ascertain the robustness of QRF for the proposed QRF-SJ model, a
comparative analysis between QRF and other ML forecasting techniques (Support
vector regression, radial basis neural network, radial basis neural network,
generalized linear model) are presented. The contributions of this paper are:
1) We implement a comprehensive probabilistic crop yield forecasting method
based on quantile random forest and Kernel density estimation (QRF-SJ). The
complete conditional probability density curve of future crop yields are
illustrated 2) Two interval prediction evaluation metrics (prediction
intervals coverage probability and prediction interval normalized average
width) are used to assess the performance of the proposed QRF-SJ. 3) We
demonstrated the superiority and feasibility of the proposed QRF-SJ model
using groundnut and millet as case studies. 4) The feature/variable importance
(a score that gives the effectiveness of each feature in predicting the yield
of the crop) is presented. This gives information to agricultural stakeholders
about the imporatant weather variable that affect the yield of the selected
crops. The rest of the paper is organized as follows: section 2 explains the
theoretical background of quantile regression, random forest, quantile random
forest, and kernel density estimation. The materials and methods are presented
in Section 3. In section 4, the results and discussion of the proposed model
using a case study of crop yield dataset is presented. The conclusion is
outlined in section 5.

##  2 Theoretical Background

This section provides a comprehensive explanation used in developing the
probabilistic crop yield forecasting. Generally, three steps are used for the
probabilistic crop yield forecasting. Firstly, the dataset is divided into a
training and testing dataset. In the second phase, the training dataset is
used to train the quantile regression forest (QRF) model. The QRF model is
then used to predict the testing data on a selected quantile. In the final
step, the probability density function are obtained by using kernel density
estimation with Epanechnikov kernel function and SJ bandwidth selection. Our
model has not been applied in other field of research or in the agriculture
sector so far. It is therefore the first of its kind in literature.

###  2.1 Quantile Regression (QR)

Conventional linear regression models make a summary of the average relation
betweeen explanatory variables $X=[X_{1},X_{2},\cdots,X_{k}]^{\prime}$ and a
response variable $Y$ depending on the conditional mean function
$\mathbb{E}(Y|X)$. It gives a partial estimate of the relationship, as it
might be needed in recounting the relationship of distinct points in the
conditional distribution of $Y$. Contrary to the conventional linear
regression, QR gives the quantiles of the conditional distribution of $Y$ as a
function of $X$ [31]. That is, QR provides much detail information about the
distribution of $Y$ than conventional linear regression model. By using QR, we
can make a good inference on the distribution of the predicted values. Machine
learning techniques that are based on quantile regression such as the quantile
random forest have an extra advantage of being able to predict non-parametric
distributions. A QR problem can be formulated as;

$q_{Y}({\tau}\mid X)=\mathbf{X}^{\prime}_{i}\bm{\beta}_{\tau}$ (1)

where $q_{Y}({\tau}\mid\cdot)$ is the conditional $\tau_{th}(0<\tau<1)$
quantile of crop yield variables $Y$, $X$ are the explanatory variables, and
$\beta_{\tau}=[\beta_{\tau}(0),\beta_{\tau}(1),\cdots\beta_{\tau}(k)]^{\prime}$
is a vector of values of quantile $\tau$. By minimizing the loss function of a
specific $\tau_{th}$ quantile, vector of values can be evaluated,

$\displaystyle\underset{\beta}{\text{min}}\sum_{i=1}^{N}\rho_{\tau}(Y_{i}-\mathbf{X}^{\prime}_{i}\bm{\beta})=$
$\displaystyle\underset{\beta}{\text{min}}\Big{[}\sum_{i:Y_{i}\geq\mathbf{X}^{\prime}_{i}\bm{\beta}}\tau\mid
Y_{i}-\mathbf{X}^{\prime}_{i}\bm{\beta}\mid+\sum_{i:Y_{i}<\mathbf{X}^{\prime}_{i}\bm{\beta}}(1-\tau)\mid
Y_{i}-\mathbf{X}^{\prime}_{i}\bm{\beta}\mid\Big{]},$ (2) $\displaystyle=$
$\displaystyle\underset{\beta}{\text{min}}\Big{[}\sum_{i}\mid\tau-\mathbf{1}_{y_{i}<\mathbf{X}^{\prime}_{i}\bm{\beta}}\mid(Y_{i}-\mathbf{X}^{\prime}_{i}\bm{\beta})\Big{]}$

Where $\mathbf{1}$ is the indicator function, $N$ is the size of the sample
data, and $\mathbf{X}_{i}=(x_{1i},x_{2i},x_{3i},\cdots,x_{ki})$ are the
independent variables. Consider the distribution of a discrete random variable
$Y_{i}$ with a less-than-well-behaved density, then the conditionbal density
function at the $\tau_{th}$ quantile given $x_{i}$ is defined as

$q_{\tau}(x)=\inf\\{y:F_{\tau}\big{(}y\mid X=x\big{)}\geq\tau\\}$ (3)

where $F_{\tau}\big{(}y\mid X=x\big{)}$ is the distribution function for
$Y_{i}$ conditional on $X_{i}$.

####  2.1.1 Prediction Interval Construction

PI is built from the conditional quantiles of the crop yield predicted from
the QRF. Particularly, the $(1-\tau)\times 100\%$ PI for crop yield ($Y$)
given the weather features ($\mathbf{X}$) is constructed as
$PI(x)=[q_{\tau/2}(Y\mid\mathbf{X}=x),q_{1-\tau/2}(Y\mid\mathbf{X}=x)]$. For
instance, the $90\%$ PI for crop yield is calculated as

$PI(x)=[q_{0.05}(Y\mid\mathbf{X}=x),q_{0.95}(Y\mid\mathbf{X}=x)].$

That is, for $x$, the crop yield is within the interval of $PI(x)$ with a very
high probability.

###  2.2 Random Forest (RF)

RF is a binary tree machine learning algorithm and a non-parametric method for
regression and classification problems. The objective of RF is to predict the
square integrable random response $Y\in\mathbb{R}$ by computing the regression
function $c(x)=\mathbb{E}\big{[}Y|{\bf X}={\bf x}\big{]}$. Assume a training
dataset $D_{n}=\\{\big{(}{\bf X}_{i},y_{i}\big{)}_{i=1}^{n}|{\bf
X}_{i}\in\mathbb{R}^{M},y\in\mathbb{R}\\}$ of an observed dataset is randomly
selected from an (unknown) probaility distribution $({\bf
x}_{i},y_{i})\sim({\bf X},Y)$. We seek to use $D_{n}$ to build an estimate.
Where ${n}$ is the total number of training samples and $M$ is the total
number of features.
Suppose $\theta$ is the parameter that determines a specific splitting node of
RF regression trees. Let $T(\theta)$ be the decision tree under consideration.
Consider the conditional distribution of $Y$ given $X=x$ depending on the
decision tree and the event that $x$ can be determined at a point on the
decision tree $R$. If there is one and only one leaf node which satisfies $x$
and is represented as $\ell(x,\theta)$ for the decision tree $T(\theta)$, then
the prediction of a single tree $T(\theta)$ for a point $x$ in the observed
data is the average over the observed values in $\ell(x,\theta)$. The weight
vector $w_{n}(x,\theta)$ for the total observation in $\ell(x,\theta)$ is
given as

$w_{n}(x,\theta)=\dfrac{\mathbf{1}_{\\{X_{n}\in
R_{\ell(x,\theta)}}\\}}{\\{p:X_{p}\in\Omega_{\ell(x,\theta)}\\}}.$ (4)

Where $\sum_{i=1}^{n}w_{n}(x,\theta)=1$, and the prediction of the single tree
$Y\mid X=x$ is the weighted average of true observation
$Y_{i}(i=1,2,\cdots,n)$,

$\hat{\vartheta}(x)=\sum_{i=1}^{n}w_{n}(x,\theta)Y_{n}$ (5)

RF uses the average prediction of $k$ individual trees, each built with an
independently and identically distributed (i.i.d.) vector
$\theta_{i}\,\,i=1,2,3,\cdots,k$ to approximate $\mathbb{E}(Y\mid X=x)$.
Denote $w_{n}(x)$ as the average of of $w_{n}(\theta)$ over the ensemble of
trees,

$w_{n}(x)=\dfrac{1}{k}\sum_{i=1}^{K}w_{n}(x,\theta_{i})$ (6)

Then, the prediction of RF is

$\hat{\vartheta}(x)=\sum_{n=1}^{N}w_{n}(x)Y_{n}$ (7)

RF estimates the conditional mean of $Y$, given $X=x$, by weighting the sum of
all the observations. The weight is larger when the conditional distribution
of $Y$ given $X=X_{n}$, is identical to the conditional distribution of $Y$
given $X=x$ [32].
RF depends on some parameters for optimal performance. The number of tress
(ntree) to grow and the number of variables that is sampled as candidates for
each split (mtry). For regression problems, mtry = $\frac{M}{3}$, where
$M$=number of features for prediction. Apart from using RF for quantile
regression forest, we shall use RF for feature importance and partial
dependence plots (PDP). For the feature importance, we use the percentage mean
decreasing accuracy (“%IncMSE") to know the importance of each of the features
in building the prediction model. The PDP illustrates how the RF model
predictions are affected by each feature assuming the rest of the features in
the RF model are controlled.

###  2.3 Quantile Random Forest (QRF)

Conventional RF predict values in individual leaf node, which is considered as
the sample mean in the leaf node. This can lead to biasness, that is extreme
values in the data samples can be over- or under-estimated. To improve the
accuracy of the prediction in the presence of extreme values in the sample
dataset, the median can be used. Hence, the median is used for point
prediction in QRF model.
QRF is a robust, non-linear, and non-parametric regression method based on
random forests method for determining conditional quantiles [29]. QRF gives an
approximation of the complete conditional distribution. Just like RF, QRF is a
set of binary regression trees. However, for each leaf node of the tree, QRF
evaluates the estimated distribution $F(y\mid X=x)=P(Y\leq y)\mid
X=x)=\mathbb{E}(1_{\\{Y\leq y\\}}\mid X=x)$ as alternative to only the mean of
$Y$ values in RF. Given a probability $p$, the quantile $q_{\tau}(X)$ is
evaluated as $\hat{q}_{\tau}(X=x_{new})=\inf\\{y:\hat{F}(y\mid
X=x_{new})>\tau\\}$. The quantiles provide a comprehensive information on the
distribution of $Y$ as a function of the predictands ($X$) than only the
conditional mean. For interval prediction,

$\Big{[}q_{\tau_{l}}(X),q_{\tau_{u}}(X)\Big{]}=\Big{[}\inf\\{y:\hat{F}(y\mid
X=x)\geq\tau_{l}\\},\inf\\{y:\hat{F}(y\mid X=x)\geq\tau_{u}\\}\Big{]}$ (8)

where $\tau_{l}<\tau_{u}$ and $\tau_{u}-\tau_{l}=\alpha$, $\alpha$ is the
probability that the predicted value fall within $y$ to lie in the interval
$[q_{\tau_{l}}(X),q_{\tau_{u}}(X)]$.
We define an approximation to the acummulated conditional probability
$\mathbb{E}(\mathbf{1}_{\\{Y\leq y\\}}\mid X=x)$ by the weighted mean of all
the observations of $\mathbf{1}_{\\{Y\leq y\\}}$ as,

$F(y\mid X=x)=\sum_{n=1}^{N}w_{n}(x)\mathbf{1}_{\\{Y_{n}\leq y\\}},$ (9)

where $w_{n}(x)$ is the same weights as in random forests. By plugging
$\hat{F}\big{(}y\mid X=x\big{)}$ into 3, the estimate $\hat{q}_{\tau}(x)$ of
the conditional quantiles $q_{\tau}(x)$ are derived,

$\hat{q}_{\tau}(x)=\hat{F}^{-1}(\tau)=\inf\\{y:\sum_{n=1}^{N}w_{n}(x)\mathbf{1}_{\\{Y_{n}\leq
y\\}}\geq p\\}$ (10)

###  2.4 Kernel density estimation using Epanechnikov Kernel function

Kernel density estimation (KDE) is a non-parametric method of estimating the
probability density function (pdf) or regression functions. KDE is basically
used for data smoothing. A Kernel density estimator at $x$ for an observed
i.i.d. and data ${\bf X}=(X_{1},X_{2},\cdots,X_{n})$ drawn from an unknown
distribution with an unknown density $f_{X}(x$), is

$\hat{f}(x;b)=\dfrac{1}{Nb}\sum_{i=1}^{N}K\Big{(}\dfrac{X_{i}-x}{b}\Big{)}$
(11)

where $N$ is the sample size, $K$ is the Kernel function, $h>0$ is the
smoothing parameter also called bandwidth. The Kernel function is non-negative
and is defined as $\int_{-\infty}^{\infty}K(x)dx=1$. Gaussian, Rectangular,
Uniform, Cosine, Epanechnikov, and Quartic are but some common examples of
kernel functions used in literatures. Different results are obtained depending
on the type of Kernel function used. In this study, we use the Epanechnikov
Kernel to build our QRF-SJ model. The Epanechnikov Kernel is defined as:

$K(u)=\begin{cases}\frac{3}{4}(1-u^{2})\qquad\mid u\mid\leq 1,\\\
0\qquad\qquad\qquad\,\,\,\text{otherwise}.\end{cases}$ (12)

Where $u=\frac{X_{i}-x}{b}$. The choice of the Epanechnikov kernel is
motivated because it has the lowest (asymptotic) mean square error (MSE) [33,
26].

#### Solve-The-Equation Plug-In Approach of Sheather and Jones (SJ)

Bandwidth determines the smoothness of the kernel density plot and is
comparable to the binwidth in a histogram. The selection of a proper bandwidth
is the most difficult problem in obtaining a good KDE [26]. A larger value
bandwidth value causes over smoothing and a very small bandwidth value causes
under smoothing. To get better results of kernel density estimator, this paper
uses Sheather and Jones (SJ) solve-the-equation (SJ) bandwidth selector [see
30, 34] for estimating the bandwidth parameter.
To quantify the accuracy of the kernel density estimator, the asymptotic mean
integrated squared error (AMISE) is used. AMISE is an approximation of mean
integrated squared error ( when $n\rightarrow\infty$ and $b=b(n)\rightarrow
0$) of $\hat{f}(x)$,

$AMISE(\hat{f}_{b}(x))=(nb)^{-1}R(K)+b^{4}R(f^{\prime\prime})\Big{(}\int
x^{2}K/2\Big{)}^{2}$ (13)

where the notation $R(g)=\int g^{2}(x)dx$ for a function $g$, $\int
x^{2}K=\int x^{2}K(x)dx$, and $f^{\prime\prime}$ is the second derivative of
$f$. The first and second term in equation 13 are the integrated variance and
integrated squared bias respectively. A very small $h$ results in a large
integrated variance and a very large $h$ results in a large integrated
squareed bias.
The flowchart of the proposed crop yield probability density forecasting
method is clearly presented in figure 1. Figure 2 shows flowchart of the
structutre of the paper.

Figure 1: The block diagram of the proposed crop yield probability density
forecasting method Figure 2: The flowchart of the structutre of the paper

##  3 Materials and Methods

###  3.1 Evaluating point prediction errors

We use root mean squared error (RMSE), mean absolute percentage error (MAPE),
R-squared ($R^{2}$), and Bias to compare the performance of different
forecasting models for point/deterministic forecasting.
RMSE (see equation 14) estimates the residual betweeen the observed and the
predicted values. The smaller the RMSE, the better the model.

$RMSE=\sqrt{\dfrac{\sum_{i=1}^{N}(\hat{y}_{i}-y_{i})^{2}}{N}}$ (14)

MAPE (see equation 15) gives the average of the absolute percentage errors. An
optimal model has the lowest MAPE.

$MAPE=\dfrac{1}{N}\sum_{i=1}^{N}\left|\dfrac{y_{i}-\hat{y}_{i}}{y_{i}}\right|\times
100\%$ (15)

$R^{2}$ (see equation 16) is a statistical measure that explains the variation
of the observed to the predicted. Generally, the closer the value of $R^{2}$
to 1 (100%), the better the model.

$R^{2}=1-\dfrac{\sum_{i=1}^{N}(y_{i}-\hat{y}_{i})^{2}}{\sum_{i=1}^{N}(y_{i}-\bar{y}_{i})^{2}}$
(16)

Bias calculates the average amount by which the observe is greater than the
predicted. For an unbiased model, the Bias should be closer to 0. Bias is
computed as in equation (17)

$Bias=\dfrac{1}{N}{\sum_{i=1}^{N}(y_{i}-\hat{y}_{i})}$ (17)

$y_{i}$, $\hat{y}_{i}$, $\bar{y}_{i}$, $N$ are the observed/actual, predicted,
mean, and the total number of dataset. The training data was used in training
the models and the testing data was used in evaluating the performance of all
the models.

###  3.2 Uncertainty of Prediction Intervals

Different metrics are used to evaluate the prediction intervals for the
results obtained from the probability density forecasting; prediction
intervals coverage probability (PICP) and prediction interval normalized
average width (PINAW).
PICP is the percentage of the testing data that fall in the interval specified
by the upper bound $U_{i}$ and the lower bound $L_{i}$ of the prediction
interval (PI). A larger PICP indicates that most of the forecasted data fall
within the PI. Generally, the value of the PICP should be greater than the
nominal confidence level.

$PICP=\dfrac{1}{N}\sum_{i=1}^{N}c_{i},$

where $N$ is the total number of years over the period of forecasting and
$c_{i}$ is a Boolean variable define as

$c_{i}=\begin{cases}1,\qquad\qquad if\,\,\,y_{i}\in[L_{i},U_{i}]\\\
0\qquad\qquad\,if\,\,\,y_{i}\notin[L_{i},U_{i}]\end{cases}$

If the quality of the forecast depends only on the PICP, the coverage
probability can be artificially improved by increasing the range between the
upper and the lower bound. However, a larger interval width is empirically not
informative. To better evaluate the quality of the PIs, the width of the PIs
must be measured. A narrow PI gives more information to the forecaster than a
wider PI. Therefore, a normalized metric PINAW which measures the average
width of the PIs can be use. PINAW is defined as:

$PINAW=\dfrac{1}{NR}\sum_{i=1}^{t}(U_{i}-L_{i})$

where $R$ is the range of the underlying targets (difference between minimum
and maximum targets)

##  4 Results and Discussion

To demonstrate the feasibility and suitability of the proposed QRF-SJ method,
historical crop yield111crop yield is defined as the harvested production of a
crop per unit of the harvested area and is measured in metric ton per hectare
(t/ha). data for different crops in Tamale metropolitan222The capital town of
the Tamale metropolitan is Tamale, which also happens to be the regional
capital of the Northern region. (a metropolis in the Northern region) were
obtained from the Statistics, Research and Information Directorate (SRID) of
the Ministry of Food and Agriculture, Ghana. The Northern region of Ghana of
which the Tamale metropolitan is located is much drier333This is because of
its close proximity to the Sahara, and the Sahel. as compared to the southern
part of Ghana and agriculture contribute more than 75% of the economic
activities in the metropolis. The region is considered to be the food basket
of Ghana. Cowpea, cassava, groundnut, maize, millet, sorghum, rice, and yam
are the major crops grown in this region. For the purpose of this research,
two of the crops (groundnut and millet) are selected as a case study to
evaluate the performance of our proposed model. These selected crops can be
used as a proxy to create an area-yield index insurance instrument for the
insurance sector. Due to the unavailability of historical data in the
metropolis, the selected crop yield data was taken from 2000 to 2016. The
growing seasons of groundnut are March-May and September-October, while the
growing season of millet is from June to October [35] . Station based daily
sunlight, humidity, precipitation, minimum temperature, maximum temperature
and average temperature from 2000 to 2016 are obtained from the Ghana
Meteorological Service. The k-nearest neighbors (KNN) algorithm was used for
imputing missing data points in the datasets. KNN locates the $k$ closest
neighbors to the observed dataset with the missing data point and imputes the
data point based on the non-missing data points in the neighbors. The average
of the climatological data over the cropping season of groundnut and millet
are computed as the growing season climatological factors.
Because of the scaling sensitivity of the inputs fed into most forecasting
techniques, the variables for the inputs are set into an identical scale. The
scale used is the min-max normalization which was set to be in the interval
$[0,1]$. The normalization is given as:

$I_{NORM}=\dfrac{I-I_{MIN}}{I_{MAX}-I_{MIN}}$

where $I_{NORM}$ is the normalized numerical value; $I_{MIN}$, $I_{MAX}$ is
the minimum and maximum values of the inputs respectively.
To validate the model, we divided the dataset into a training (80%) and
testing (20%) dataset. That is, the dataset from 2000-2013 are selected as the
training data and 2014-2016 are selected as the testing data.

$R$ statistical software was used in implementing all the forecasting
approaches.

Crop | Mean | Std | Min | Max | Skewness
---|---|---|---|---|---
Groundnut yield (t/ha) | 1.22 | 0.50 | 0.50 | 1.90 | 0.05
Millet yield (t/ha) | 1.19 | 0.30 | 0.72 | 1.70 | 0.10
Table 1: Summary Statistics of Groundnut and Millet Yield

###  4.1 Empirical results and analysis of the Models

In order to prove the superiority of the QRF, it is compared with some popular
forecasting techniques like Radial Basis Neural Network (NN), Generalized
Linear Model (GLM), Support Vector Regression with linear (SVR (linear)) and
radial basis (SVR (radial)) kernel function. The hyperparameters used in
tuning the ML techniques is detailed in appendices LABEL:appendix_groundnut
and LABEL:appendix_millet.
The evaluation metrics (RMSE, MAPE, $R^{2}$, and Bias) of these methods are
given in Table 2. Figures 3 and 4 show the visual performance of the
evaluation metrics of the forecasting techniques. Comparative to the four
benchmark methods, it is evident that QRF (both mean and median prediction)
performed better in terms of RMSE, MAPE and Bias when predicting the yield of
groundnut and millet for the testing data. GLM outperformed all the other
methods when using $R^{2}$ as the performance metrics for both groundnut and
millet. It can be observed that the QRF median prediction is the same as the
mean prediction. Quantitatively, the RMSE, MAPE, $R^{2}$ and Bias of QRF (mean
and median) for groundnut are $0.3787$ t/ha, $24.0026\%$, $0.9830$, and
$0.3394$ t/ha respectively. The RMSE, MAPE, $R^{2}$ and Bias of QRF (mean and
median) for millet are $0.0173$ t/ha, $0.9090\%$, $0.9805$ and $0.01$ t/ha
respectively. Generally, QRF is optimal in predicting the yield of all the two
crops. For this reason, we conclude that QRF is the best forecasting technique
for predicting the yield of the two crops as compared to the other benchmark
forecasting techniques. This motivated us to use QRF for the probability
density forecasting.

Method | Groundnut |  Millet
---|---|---
| RMSE (t/ha) | MAPE (%) | $R^{2}$ | Bias (t/ha) |  | RMSE (t/ha) | MAPE (%) | $R^{2}$ | Bias(t/ha)
QRF(mean) |  | 0.3787 | 24.0026 | 0.9830 | 0.3394 |  | 0.0173 | 0.9090 | 0.9805 | 0.0100
QRF(median) |  | 0.3787 | 24.0026 | 0.9830 | 0.3394 |  | 0.0173 | 0.9090 | 0.9805 | 0.0100
SVR (radial) |  | 0.4190 | 29.8431 | 0.7050 | 0.4133 |  | 0.1035 | 7.8520 | 0.0166 | 0.0046
SVR (linear) |  | 0.7262 | 52.5940 | 0.5959 | 0.7198 |  | 0.1559 | 13.9280 | 0.2407 | 0.0691
NN |  | 0.6707 | 48.4656 | 0.5575 | 0.6679 |  | 0.2860 | 25.3258 | 0.8379 | 0.2827
GLM |  | 0.5043 | 36.7089 | 0.9962 | 0.4989 |  | 0.3033 | 27.2320 | 0.9983 | 0.2945

Table 2: Evaluation metrics (RMSE, MAPE, R-squared, Bias) of point prediction
via QRFs, SVR(radial), SVR(linear), NN, and GLM using testing data.

Figure 3: Bar chart of RMSE, MAPE, R-Squared, and Bias of observed and
predicted groundnut yield by QRFs, SVR(radial), SVR(linear), NN, and GLM using
testing data. QRF (mean): mean prediction of crop yield given weather features
using QRF; QRF (median): median prediction of crop yield response using QRF;
SVR(radial): prediction of crop yield using SVR radial basis kernel function;
SVR(linear): prediction of crop yield using SVR linear kernel.

Figure 4: Bar chart of RMSE, MAPE, R-Squared, and Bias of observed and
predicted millet yield by QRFs, SVR(radial), SVR(linear), NN, and GLM using
testing data. QRF (mean): mean prediction of crop yield given weather features
using QRF; QRF (median): median prediction of crop yield response using QRF;
SVR(radial): prediction of crop yield using SVR radial basis kernel function;
SVR(linear): prediction of crop yield using SVR linear kernel.

Table 3 presents the values of the prediction results and prediction intervals
based on QRF-SJ forecast for testing data. From this table, it is clear that
each prediction has its own uncertainty bound. Figure 5 gives the point
forecast and the prediction interval of groundnut and millet respectively. To
assess the performance of point forecast for the model, the performance
metrics over the predicted period for the two crops are presented in table 4.
It is clear from figure 5 that the target values lie within the prediction
interval. We can therefore conclude that the proposed model captures the
uncertainty of the two crop yields accurately.

Year | Groundnut |  Millet
---|---|---
| Lower bound | Observed | Predicted | Upper bound |  | Lower bound | Observed | Predicted | Upper bound
2014 |  | 0.5630 | 1.5000 | 1.1100 | 1.9500 |  | 0.7420 | 1.2400 | 1.2400 | 1.6101
2015 |  | 0.5667 | 1.3300 | 1.0700 | 1.8800 |  | 0.7397 | 1.0000 | 1.0000 | 1.6007
2016 |  | 0.5680 | 1.3000 | 1.1200 | 1.9200 |  | 0.7429 | 1.1000 | 1.0700 | 1.6120

Table 3: Prediction results and prediction intervals for the yield of
groundnut and millet for testing data

Crop | Confidence level (%) | MAPE (%) | RMSE (t/ha) | Bias (t/ha)
---|---|---|---|---
Groundnut | 90 | 27.50 | 0.2895 | 0.2399
Millet | 90 | $1\times 10^{-14}$ | 0.0072 | 0.0147

Table 4: Evaluation metrics of point prediction via QRF-SJ for testing data

Figure 5: Prediction results and prediction intervals for the yield of
groundnut and millet.

To show the satisfactory performance of the proposed QRF-SJ model, PICP and
PINAW are used as the evaluation metrics. The performance of the model is
presented in Table 5. It is clear from the table that, all observed values
were captured by the confidence interval. This is very important for effective
decision-making process. For both groundnut and millet, the PICP was evaluated
to be 100%. The PINAW of groundnut is however smaller than that of millet.
Considering the high variance of the yield of groundnut in table 1, the
probabilistic performance of our proposed method is sufficient.
A visual representation of the probability density curve for the predictions
based on QRF-SJ for the yield of groundnut and millet is presented in figures
6 and 7 respectively. The actual crop yield for the specific year are
presented in orchid, blue, and red dots for 2014, 2015, and 2016 respectively.
In figures 6 and 7, the respective actual crop yields of groundnut and millet
for each of the predicted year fall within the predicted region of the
forecast distribution. The probability density curve gives the complete
probability distribution of the future crop yield and hence the uncertainty of
the forecasting can be quantified.

Crop | Confidence level (%) | PICP (%) | PINAW (%)
---|---|---|---
Groundnut | 90 | 100 | 12.65
Millet | 90 | 100 | 16.76
Table 5: Prediction Interval evaluation indices of QRF-SJ method Figure 6:
probability density curve based on QRF-SJ from 2014-2016, the dots on the
x-axis represents the actual values of the crop yield. Figure 7: probability
density curve based on QRF-SJ from 2014-2016, the dots on the x-axis
represents the actual values of the crop yield.

###  4.2 Feature Importance

The level of feature importance measures according to random forest is shown
in Table 6. The higher the percentage increase in mean square error (%IncMSE)
of a feature, the higher the importance of that feature in the prediction
model.
From the feature importance measure in Table 6, average temperature, minimum
temperature, and rainfall are the three most important features among the six
features that affect the yield of groundnut. Humidity, rainfall, and average
temperature are the top three features that influences the yield of millet.
The amount of sunshine does not have a lot of effect on the yield of
groundnut. Maximum temperature do not have a lot of effect on the yield of
both crops. Generally, average temperature do have a lot of effect on the
yield of these two crops. The partial dependence plots (PDP) in Figures 8 and
10 show the marginal effect of the top three important features of the yield
of groundnut and millet respectively.
From the PDP of groundnut, an increase in the amount of average and minimum
temperatures result in a decrease in the yield of groundnut. However, an
increase in the amount of rainfall increases the yield of groundnut. An
increasing relative humidity and minimum temperatures decreases the yield of
millet. In both crops, the yield is minimum when the minimum temperature is
around 23.2∘C. From the PDPs, it can be concluded that the yield of groundnut
and millet increases as the amount of rainfall increases to about 100mm. The
yield of groundnut and millet decreases as the minimum temperature increases
from about 22∘C to 23.2∘C. Figures 9 and 11 show the three-dimensional (3-D)
partial dependence of the top 3 ranked features from feature importance
measures of Random Forests model for groundnut and millet respectively. From
table 6, it is clear that all the six weather features do influence the yield
of all the two crops.

Feature | Groundnut |  Millet
---|---|---
Rank | %IncMSE |  | Rank | %IncMSe
Sunshine | 6 | 0.101 |  | 4 | 0.120
Humidity | 4 | 0.144 |  | 1 | 0.210
Rainfall | 3 | 0.162 |  | 2 | 0.198
AvgT | 1 | 0.175 |  | 3 | 0.121
MaxT | 5 | 0.125 |  | 5 | 0.114
MinT | 2 | 0.171 |  | 6 | 0.115
Table 6: Rank corresponds to variable importance measure determined by Random
Forest (RF) model for each crop dataset. AvgT, MaxT, MinT represent average,
maximum and minimum temperature respectively. Figure 8: Partial dependence
plot of groundnut for the top 3 ranked predictor variable from variable
importance measures of Random Forests models using the testing data Figure 9:
Groundnut: A 3-D surface Partial dependence of Rainfall on Average and Minimum
Temperature based on a random forest using the testing data. Figure 10:
Partial dependence plot of millet for the top 3 ranked predictor variable from
variable importance measures of Random Forests models using the testing data.
Figure 11: Millet: A 3-D surface Partial dependence of Average temperature on
humidity and rainfall based on a random forest using the testing data.

###  4.3 Discussion

####  4.3.1 QRF-SJ method

Due to the variability of the yield of crops as a result of climate
variability [36], a method QRF-SJ that can effectively handle this effect is
proposed using quantile random forest and kernel density estimation. QRF does
not only outperform other benchmark methods like SVR, NN, and GLM in point
forecasting but also gives useful probability distribution when combined with
kernel density estimation. Probability density forecasting has the ability to
represent uncertainty of crop yields as the probability distributions around a
prediction interval. Forecasting the distributions serves as an indicator for
the accuracy of forecast and provides important information for decision
making.
Crop yield forecasting is an integral factor in precision farming and can
promote the expansion of the agriculture sector. On the part of the
government, this forecasting model will aid in effective decision planning.
Hence, food shortages can be avoided and if possible governments can arrange
for food imports rather than seeking for emergency food assistance. The
proposed forecasting model can be used for trade development policies and
other humanitarian assistance connected to food security. For the industry
players like the insurance and financial sector, crop yield forecasting helps
in measuring crop loss’ in advance. Consequently, fair premium rates and
pricing of agricultural index insurance and weather derivatives can be
determined. The farmer however is able to measure the future uncertainties of
the farm produce and make effective plans for a set of possible outcomes and
in particularly precision farming.

####  4.3.2 Feature Importance

The feature importance gave a score of the importance of each weather variable
in building the QRF model. By identifying the contributions of these weather
variables to crop yields and applying effective forecasting models, the farmer
for instance will be able to detect the specific weather variable that affects
the yield of the crops and can adapt the appropriate risk management
strategies. The insurance sector is able to realize the specific weather
variable which correlates with the actual farm yield. Using the feature
importance, insurance companies can sell different weather-index insurance
products to farmers depending on the weather variable that affect the yield of
their crops. Payouts of this weather insurance will be triggered based on the
level for which the weather variable affects the yield of the crop. Base on
feature importance, insurers and re-insurers companies are able to price their
insurance product base on the specific climatic factors that affect the yield
of a specific crop. They can also merge the most important climatic factors
(compouns events) in their pricing models.

##  5 Conclusion

Crop yield forecasting that considers the uncertainty of weather as a result
of changes in climate is of crucial importance for efficient decision-making
process in the agricultural sector. In this paper, a hybrid crop yield
probability density forecasting method that has the potential to draw complete
conditional probability density curve of future crop yields is presented. In
the method, quantile regression forest is used to build the nonlinear quantile
regression forecasting model and to capture the nonlinear relationship between
the weather variables and crop yields. Epanechnikov kernel function and solve-
the equation plug-in approach of Sheather and Jones are employed in the method
to construct the probability density forecasting curve. Prediction interval
coverage probability and prediction interval normalized average width are used
to evaluate the quality of the prediction intervals constructed by QRF-SJ. The
performance and accuracy of the QRF-SJ crop yield forecasting model are
evaluated using two real dataset (groundnut and millet yields) as case
studies. The numerical results gave a 100% PICPs and narrow PINAWs. All the
observed groundnut and millet crop yield values are located in the probability
density curves. The results show the superiority and feasibility of the
proposed QRF-SJ model in forecasting the yield of crops in the midst of of
weather uncertainties. Using the feature importance, agricultural stakeholders
will be able to detect the major weather variables that affect the yield of
crops and can make pragmatic interventions to curtail any uncertainties.

## Acknowledgements

The first author wishes to thank African Union and Pan African University,
Institute for Basic Sciences Technology and Innovation, Kenya, for their
financial support for this research.

## References

  * Rana and Rana [2014] S. S. Rana, R. S. Rana, ADVANCES IN CROP GROWTH AND PRODUCTIVITY, Technical Report, Publication of the Department of Agronomy, CSK Himachal Pradesh Krishi Vishvavidyalaya, Palampur, India, 2014\.
  * Brevik [2013] E. Brevik,  The potential impact of climate change on soil properties and processes and corresponding influence on food security,  Agriculture 3 (2013) 398–417.
  * Chen et al. [2004] C.-C. Chen, B. A. McCarl, D. E. Schimmelpfennig,  Yield variability as influenced by climate: A statistical investigation,  Climatic Change 66 (2004) 239–261.
  * Isik and Devadoss [2006] M. Isik, S. Devadoss,  An analysis of the impact of climate change on crop yields and yield variability,  Applied Economics 38 (2006) 835–844.
  * Southworth et al. [2000] J. Southworth, J. Randolph, M. Habeck, O. Doering, R. Pfeifer, D. G. Rao, J. Johnston,  Consequences of future climate change and changing climate variability on maize yields in the midwestern united states,  Agriculture, Ecosystems & Environment 82 (2000) 139–158.
  * Asfaw and Lipper [2016] S. Asfaw, L. Lipper, Managing climate risk using climate-smart agriculture, FAO, 2016.
  * Gornall et al. [2010] J. Gornall, R. Betts, E. Burke, R. Clark, J. Camp, K. Willett, A. Wiltshire,  Implications of climate change for agricultural productivity in the early twenty-first century,  Philosophical Transactions of the Royal Society B: Biological Sciences 365 (2010) 2973–2989.
  * Gyamerah et al. [2018] S. A. Gyamerah, P. Ngare, D. Ikpe,  Regime-switching temperature dynamics model for weather derivatives,  International Journal of Stochastic Analysis 2018 (2018).
  * Delincé et al. [2017] J. Delincé, et al.,  Recent practices and advances for amis crop yield forecasting at farm and parcel level: a review.,  Recent practices and advances for AMIS crop yield forecasting at farm and parcel level: a review. (2017).
  * Shi et al. [2013] W. Shi, F. Tao, Z. Zhang,  A review on statistical models for identifying climate contributions to crop yields,  Journal of geographical sciences 23 (2013) 567–576.
  * Michler et al. [2015] J. D. Michler, F. G. Viens, G. E. Shively, et al.,  Risk, agricultural production, and weather index insurance in village south asia,  in: International Association of Agricultural Economists 2015 International Conference of Agricultural Economists, Milan, Italy, August, pp. 8–14.
  * Lobell and Burke [2010] D. B. Lobell, M. B. Burke,  On the use of statistical models to predict crop yield responses to climate change,  Agricultural and Forest Meteorology 150 (2010) 1443–1452.
  * Choudhury and Jones [2014] A. Choudhury, J. Jones,  Crop yield prediction using time series models,  Journal of Economics and Economic Education Research 15 (2014) 53.
  * Hong-ying et al. [2008] L. Hong-ying, H. Yan-lin, Z. Yong-juan, et al.,  Variations trend of grain yield per unit area and fertilizer application systems in liaoning province,  System Sciences and Comprehensive Studies in Agriculture 24 (2008) 408–410.
  * ZHANG and ZHANG [2007] X.-j. ZHANG, X.-l. ZHANG,  Application of time series analysis model on total corn yield of shandong province [j],  Research of Soil and Water Conservation 3 (2007) 098.
  * Yu and Babcock [2011] T. Yu, B. A. Babcock,  Estimating non-linear weather impacts on corn yield—a bayesian approach (2011).
  * Chmielewski and Potts [1995] F.-M. Chmielewski, J. M. Potts,  The relationship between crop yields from an experiment in southern england and long-term climate variations,  Agricultural and Forest Meteorology 73 (1995) 43–66.
  * Jeong et al. [2016] J. H. Jeong, J. P. Resop, N. D. Mueller, D. H. Fleisher, K. Yun, E. E. Butler, D. J. Timlin, K.-M. Shim, J. S. Gerber, V. R. Reddy, et al.,  Random forests for global and regional crop yield predictions,  PLoS One 11 (2016) e0156571.
  * Everingham et al. [2016] Y. Everingham, J. Sexton, D. Skocaj, G. Inman-Bamber,  Accurate prediction of sugarcane yield using a random forest algorithm,  Agronomy for sustainable development 36 (2016) 27.
  * Dahikar and Rode [2014] S. S. Dahikar, S. V. Rode,  Agricultural crop yield prediction using artificial neural network approach,  International Journal of Innovative Research in Electrical, Electronics, Instrumentation and Control Engineering 2 (2014) 683–686.
  * Fang et al. [2018] Y. Fang, P. Xu, J. Yang, Y. Qin,  A quantile regression forest based method to predict drug response and assess prediction reliability,  PloS one 13 (2018) e0205155.
  * Jiang et al. [2017] F. Jiang, W. Wu, Z. Peng,  A semi-parametric quantile regression random forest approach for evaluating muti-period value at risk,  in: Control Conference (CCC), 2017 36th Chinese, IEEE, pp. 5642–5646.
  * Chatfield [2000] C. Chatfield, Time-series forecasting, Chapman and Hall/CRC, 2000.
  * Deser et al. [2012] C. Deser, A. Phillips, V. Bourdette, H. Teng,  Uncertainty in climate change projections: the role of internal variability,  Climate dynamics 38 (2012) 527–546.
  * Reilly et al. [2001] J. Reilly, P. H. Stone, C. E. Forest, M. D. Webster, H. D. Jacoby, R. G. Prinn, Uncertainty and climate change assessments, 2001\.
  * Wand and Jones [1994] M. P. Wand, M. C. Jones, Kernel smoothing, Chapman and Hall/CRC, 1994\.
  * Barnwal and Kotani [2013] P. Barnwal, K. Kotani,  Climatic impacts across agricultural crop yield distributions: An application of quantile regression on rice crops in andhra pradesh, india,  Ecological Economics 87 (2013) 95–109.
  * Wang [2012] J. Wang,  Bayesian quantile regression for parametric nonlinear mixed effects models,  Statistical Methods & Applications 21 (2012) 279–295.
  * Meinshausen [2006] N. Meinshausen,  Quantile regression forests,  Journal of Machine Learning Research 7 (2006) 983–999.
  * Jones et al. [1996] M. C. Jones, J. S. Marron, S. J. Sheather,  A brief survey of bandwidth selection for density estimation,  Journal of the American statistical association 91 (1996) 401–407.
  * Koenker and Bassett [1978] R. Koenker, G. Bassett,  Regression quantiles. econometrica 46 33–50,  Mathematical Reviews (MathSciNet): MR474644 Digital Object Identifier: doi 10 (1978) 1913643\.
  * Lin and Jeon [2006] Y. Lin, Y. Jeon,  Random forests and adaptive nearest neighbors,  Journal of the American Statistical Association 101 (2006) 578–590.
  * Epanechnikov [1969] V. A. Epanechnikov,  Non-parametric estimation of a multivariate probability density,  Theory of Probability & Its Applications 14 (1969) 153–158.
  * Sheather and Jones [1991] S. J. Sheather, M. C. Jones,  A reliable data-based bandwidth selection method for kernel density estimation,  Journal of the Royal Statistical Society: Series B (Methodological) 53 (1991) 683–690.
  * FAO [2010] FAO, Crop calendar - an information tool for seed security, $http://www.fao.org/agriculture/seed/cropcalendar/welcome.do$, 2010\. (accessed 26 September 2018).
  * Alexandrov and Hoogenboom [2000] V. Alexandrov, G. Hoogenboom,  The impact of climate variability and change on crop yield in bulgaria,  Agricultural and forest meteorology 104 (2000) 315–327.