appendixfiguresapfList of Figures in Supplementary Material
appendixtablesaptList of Tables in Supplementary Material

# Human limits in Machine Learning: Prediction of plant phenotypes using soil
microbiome data

Rosa Aghdam
Wisconsin Institute for Discovery
University of Wisconsin-Madison
Madison, WI

&Xudong Tang∗
Department of Statistics
Wisconsin Institute for Discovery
University of Wisconsin-Madison
Madison, WI

&Shan Shan
Department of Plant Pathology
University of Wisconsin-Madison
Madison, WI
&Richard Lankau
Department of Plant Pathology
University of Wisconsin-Madison
Madison, WI
&Claudia Solís-Lemus
Department of Plant Pathology
Wisconsin Institute for Discovery
University of Wisconsin-Madison
Madison, WI
Joint first authors with equal contribution in alphabetical orderCorresponding
author<EMAIL_ADDRESS>
###### Abstract

The preservation of soil health has been identified as one of the main
challenges of the XXI century given its vast (and potentially threatening)
ramifications in agriculture, human health and biodiversity. Here, we provide
the first deep investigation of the predictive potential of machine-learning
models to understand the connections between soil and biological phenotypes.
Indeed, we investigate an integrative framework performing accurate machine-
learning-based prediction of plant phenotypes from biological, chemical and
physical properties of the soil via two models: random forest and Bayesian
neural network. We show that prediction is improved, as evidenced by higher
weighted F1 scores, when incorporating into the models environmental features
like soil physicochemical properties and microbial population density in
addition to the microbiome information. Furthermore, by exploring multiple
data preprocessing strategies such as normalization, zero replacement, and
data augmentation, we confirm that human decisions have a huge impact on the
predictive performance. In particular, we show that the naive total sum
scaling normalization that is commonly used in microbiome research is not the
optimal strategy to maximize predictive power. In addition, we find that
accurately defined labels are more important than normalization, taxonomic
level or model characteristics. That is, if humans are unable to classify the
samples and provide accurate labels, the performance of machine-learning
models will be limited. Lastly, we present strategies for domain scientists
via a full model selection decision tree to identify the human choices that
maximize the prediction power of the models. Our work is accompanied by open
source reproducible scripts (https://github.com/solislemuslab/soil-microbiome-
nn) for maximum outreach among the microbiome research community.

##  Appendix A Introduction

It is undeniable that machine learning (ML) has revolutionized the manner in
which scientific research is performed in recent years. Among the many tasks
performed by ML models in our daily lives, researchers have relied on ML to
assist in clinical diagnoses [13], identification of bacterial phenotypes such
as antimicrobial resistance [3], and even identification of objects in space
[15]. Recently, the vast evidence of the prediction power of ML models on a
wide range of applications has launched the adoption of these models in other
domains such as sustainable agriculture where soil health – characterized by a
wide range of biological, chemical, and physical properties [28] – is explored
as the potential driver to predict plant and crop phenotypes, such as disease
susceptibility or yield.

The question of whether the improvement of soil health could result in
superior crop yield and disease resistance remains open, as researchers have
not been able to identify a set of indicators to accurately and robustly
predict plant outcomes from soil information. Among all the candidate
indicators for soil health, soil microbiome is one that continues to be
understudied in its predictive potential in the productivity and resilience of
agricultural ecosystems [27]. Aided by amplicon sequencing of highly preserved
phylogenetically informative marker genes, like the 16S ribosomal RNA for
bacteria and the internal transcribed spacer (ITS) for fungi, there have been
extensive studies on the complexity and diversity of soil microbial
communities over the last decade, yet little is still understood on how
changes in the microbial community in soil would directly impact on plant
growth and health. While it is already well-recognized that soil microbiome is
closely related to plant health and productivity [6], most current research
efforts on the soil microbiome are limited to the change on microbial
community in response to agricultural management.

Fortunately, the rapid development of ML models could fill this gap. Over the
years, ML models have been used on high-dimensional, highly interconnected
data for prediction and decision-making that does not rely on a priori
knowledge on the effects and interactions among covariates. Since we have
little or no knowledge of most species among the thousands contained in the
soil microbial communities, it would be beneficial to exploit the power of ML
models on microbiome data as they are able to explore the unknown interactions
between the microbial communities and crop phenotypes.

There are three main challenges, however, for analyzing microbiome datasets
with ML methodologies: (1) Data are inherently compositional, which means that
the raw counts are normalized based on the total number of reads that were
collected from the sample. Therefore, microbial abundances are not
independent, and the use of traditional statistical techniques (such as
correlation) might lead to increased false discovery rates [11]. (2) Data are
highly sparse, which means that datasets include a large number of operational
taxonomic units (OTUs) that are present in a small proportion of samples (or
in none at all) [14]. (3) Data are high-dimensional, which means that the
number of OTUs are larger than the number of samples, especially in more
specific taxonomic orders like Order, Family, and Genus.

In addition to challenges on the microbiome data, prediction power also
deteriorates on settings involving inaccurate and imbalanced labels. That is,
when we consider the task of prediction, each microbial sample should be
labeled by a plant outcome value (say, high yield or low yield). It is common
that labels in these datasets are imbalanced with one class representing the
majority of observed samples which is denoted as an imbalanced binary
classification problem. Furthermore, decisions on the class labels are many
times not straight-forward. While diseased plant or non-diseased plant tends
to be an indisputable classification, determining what constitutes high versus
low yield is up for debate and depends not only on the season and other
environmental factors, but also on the specific crop variety.

Here, we explore the potential of ML models in the prediction of crop
phenotypes of yield and disease from soil data. We utilize a real dataset from
potato fields in Wisconsin and Minnesota and focus on the performance of the
models while facing data challenges related to binarization, imbalance,
compositionality, sparsity and high dimensionality. Furthermore, we test the
impact of specific data preprocessing steps such as 1) different normalization
and zero replacement strategies to overcome the compositionality and sparsity,
2) different feature selection strategies to overcome the high dimensionality,
and 3) data augmentation techniques to overcome the imbalance of the data. In
addition, we are also interested in answering the question of whether soil
microbiome data alone has enough predictive power to predict the disease and
yield phenotypes or whether other information on the soil, such as chemical
composition, is necessary for accurate prediction. The answer to this question
will inform farmers on the best data collection strategies given that soil
microbiome data is expensive. That is, if the soil microbiome data does not
provide enough predictive power compared to other (cheaper) soil measurements,
then the collection of microbiome data for prediction purposes would be
futile.

For our investigation, we select two distinct ML models: 1) Random forests
(RF) as they have always been identified as the top predicting models in a
variety of domain applications [5]. 2) Neural networks given that they are
capable of exploring the arbitrary complex interactions among features that we
do not have prior knowledge of. However, with the very small sample size and
large feature space of this study, the traditional back-propagation-based
neural network models will likely have a large bias or overfit [12]. Thus, we
select the Bayesian neural network (BNN) as our second model, which offers
natural protection against overfitting [9].

Among the main findings, we can highlight that microbiome data alone indeed
has predictive power for disease outcomes, especially for pitted scab disease,
but not to predict yield. We also find that the most powerful prediction is
achieved by combination of environmental information and microbiome data.
Among the data preprocessing strategies that we explored, we find that
normalization and zero replacement strategies have a huge impact on the
prediction power of the models, yet there are strong interaction effects with
taxonomic levels, and thus, it is impossible to identify one strategy that
outperforms the others. In terms of the model, we found no clear differences
between the RF and the BNN, yet the latter is more computationally expensive
and prohibited on some datasets with a large number of predictors, which
prompts us to prefer RF for most settings. We conclude our investigation with
a full model selection (FMS) decision tree [25, 26] that allows us to identify
combinations of normalization, zero replacement, feature selection, and model
that yield the highest prediction accuracy and that can serve to provide
specific strategies for other researchers with similar datasets.

##  Appendix B Results

Figure 1 shows a graphical representation of our pipeline with three major
steps: (1) data preparation, (2) feature selection, and (3) classification
based on random forest (RF) and Bayesian neural network (BNN) models with
various types of predictors, including microbiome data, environmental, and a
combination of both. Every panel in this graphical abstract will be described
in more details in the Methods section.

Figure 1: Workflow of the analyses with three main steps: (1) Data
Preparation, (2) Feature Selection, and (3) Classification. In (1) Data
Preparation, we consider OTUs at different taxonomic levels and filter by
sample size. In addition, we perform five normalization methods and four zero
replacement methods. In (2) Feature selection, we rank OTUs based on the
number of times they are selected as important features by machine learning
(ML) criteria. Then, we select OTUs that display the greatest degree of
difference on microbial networks reconstructed from samples of each class.
Finally, we score OTUs based on whether they are selected by ML ($score=1$),
by network comparison ($score=2$), both ($score=3$) or neither ($score=0$). In
(3) Classification, the Venn diagram depicts the different types of
predictors: microbiome (OTUs), environmental (Env), and the combination of
both. random forest and Bayesian neural network classification models are
fitted on the different input predictors.

###  B.1 Overall performance of predictive models: Manual binarization causes
inaccurate prediction of yield

First, we identify the outcomes (disease or yield) that are accurately
predicted across models (and are thus robust to prediction regardless of model
choices), as well as the models that accurately predict across outcomes (and
are thus the most powerful model alternatives). To do so, we aggregate the
weighted F1 scores on data preprocessing choices such as normalization, zero
replacement, and taxonomic levels for every model and every outcome. Figure S3
(random forest) and Figure S4 (BNN) show the boxplots with the weighted F1
scores of the 14 models described in Table 3. Columns correspond to the six
responses: four diseases and two yield outcomes. For a given panel (model in
row and response in column), the boxplot corresponds to the different weighted
F1 scores for every combination of normalizations/zero replacement strategies
as well as different taxonomic levels (Table 4). For example, the boxplots for
the ALL-OTU model (first row) include weighted F1 scores of the model fit on
20 normalization/zero replacement strategies, and 5 taxonomic levels (100
different weighted F1 scores per outcome). The dashed line in each panel
corresponds to the average weighted F1 score of the model when fit with all
random datasets (see Section B.2.1). This line allows us to assess whether the
real data has more predictive power than random data.

While these plots do not allow us to distinguish differences by taxonomic
level or normalization/zero replacement strategy (more on that in the next
subsections), we can identify outcomes (columns) that can be more accurately
predicted across models (rows). Additionally, we can identify models that are
capable of accurately predicting more outcomes. It is readily evident, for
example, that the yield outcomes cannot be accurately predicted by any model
as all the weighted F1 scores fall consistently below dashed line. It is
notable that the yield responses are those for which there is not a clear
binarization strategy. Since we are artificially separating samples into the
two classes (low and high yield) based on whether they are above or below the
variety-specific median, samples on the boundary will in fact be very similar
to each other, and thus, difficult to classify. Furthermore, the poor
prediction of yield is not restricted to one data type (microbiome vs
environmental) which also suggests that the prediction challenges arise from
the binarization process rather than the model or set of predictors.

Disease outcomes, on the contrary, display higher weighted F1 scores overall,
and in particular, pitted scab displays weighted F1 scores that are
consistently above the random prediction dashed line across different models.
For the case of black scurf, even when the weighted F1 scores are consistently
above dashed line, this is a deceiving result, as this disease outcome is
highly imbalanced. This means that a naive model predicting all samples to
belong to the majority class will have high prediction accuracy (see dashed
line above 0.8 for random data). We investigate the prediction of black scurf
more carefully with the augmented data that balances the proportion of both
classes (Section B.2.4). Finally, given the poor performance on yield, we
focus on the fine-grained description of results for the disease outcomes only
for the remaining of the manuscript.

Regarding differences by model (rows), the most accurate prediction across all
outcomes can be obtained with alpha diversity and soil chemistry data
(Alpha+Soil) for the case of RF (Figure S3) and with OTUs identified as
important by both ML and network comparison strategies (more on feature
selection in Section B.2.3) along with soil chemistry data (OTU-S3+Soil) for
the case of BNN models (Figure S4). Both of these results show that
incorporating environmental information in conjunction with microbiome data
can increase the predictive power for diseases outcomes.

###  B.2 Prediction of potato disease from microbiome data

####  B.2.1 More power to predict pitted scab disease in microbiome data
compared to random data

In order to study whether the microbiome data has predictive signal to
classify samples corresponding to diseased or non-diseased potatoes, we can
compare the performance of our models when fitted on the real data and when
fitted on completely random data. Intuitively, if the data has predictive
signal, models fitted on real data will dramatically outperform models fitted
on random data.

Figure 2 shows the scatterplot of F1 scores in the samples in one class (non-
diseased potato: Label $0$) versus the F1 scores in the samples in the other
class (diseased potato: Label $1$) for two disease outcomes (rows) and the
four strategies to generate random data (columns). Strategies correspond to
randomization of rows or columns, or generation of random abundance values
(more details on the randomization strategies in Methods). Each blue point in
this scatterplot corresponds a random dataset, and the red point corresponds
to the F1 scores when the model is fitted with the real data. Here, we only
focus on microbiome data at the Phylum level.

For the pitted scab on the top row, the F1 scores for the original dataset
(red dots) are on the upper right quadrant compared to points from random
datasets (blue dots). This means that the F1 scores on the original data are
higher than F1 scores on random data, and this, the original data has more
predictive power than random data. However, for the superficial scab on the
bottom row, many random datasets show better performance (higher F1 score
values) than the real data, and thus, the real data does not contain much
information to predict this disease. Other diseases as well as yield outcomes
also underperform compared to random data as shown in Figure S5 in the
Appendix. We highlight that pitted scab is the most severe case of the potato
disease common scab, and possibly a stronger signal of the soil microbial
community composition and corresponding disease-suppressiveness than other
measures of the disease such as superficial scab.

Figure 2: Scatterplot for F1 scores in the non-diseased (Label $0$) group
versus the F1 scores in the diseased (Label $1$) group by disease outcome
(rows) and randomization strategy (columns). Each panel has 200 generated
datasets (each blue point is a random dataset). The red point corresponds to
the F1 scores when the model is fitted on the real data. These results are for
OTUs at the Phylum level. There is clear predictive power for the pitted scab
(Scabpit) as the red dot is in the upper right region of the panel, but none
for the superficial scab (Scabsuper).

####  B.2.2 No normalization or zero replacement strategy proven Superior in
terms of predictive power

One of the goals of our study is to identify ideal data preprocessing steps
that are guaranteed to maximize predictive power on the machine learning
models. Our results, however, suggest that there is not a clear pattern of
superiority among the normalization or zero replacement strategies, and that
selection of the data preprocessing steps needs to be data specific (and
taxonomic level specific).

Figure 3 shows the weighted F1 scores for pitted scab for different
combinations of normalization and zero replacement strategies (x-axis) for the
two types of models (RF and BNN). These analyses include all OTUs under the
five taxonomic levels (different colors). Similar plots for other diseases are
presented in Figures 6 to S10 in the Appendix. As mentioned, we note that
there is not an apparent pattern that points at any normalization or zero
replacement strategy being superior to other alternatives across all taxonomic
levels. In fact, there is considerable interaction between the
normalization/zero replacement method and the taxonomic level. For example,
for the RF model, the best result is achieved with Phylum level and cumulative
sum scaling normalization with pseudo-zero replacement strategy (CSS+pseudo)
or common sum scaling normalization without any zero replacement strategy
(COM+none). For BNN, however, the best results are achieved with the common
sum scaling normalization with the multiplicative zero replacement
(COM+multRepl) for the Phylum level.

For a given normalization/zero replacement strategy (x-axis), the fact that
the points in the scatterplot are not overlapping shows that taxonomic levels
have a impact on the predictive power of the model. There is, however, no
apparent pattern on a given taxonomic level providing more predictive signal.
Phylum level seems to produce higher weighted F1 scores for some
normalization/zero replacement strategies, but not all.

When we compare the range of weighted F1 scores across normalization and zero
replacement strategies, we see that the effect of the strategy is not
negligible. For example, for a Phylum level, the lowest weighted F1 score is
around 0.75 for centered log-ratio normalization with pseudo-zero replacement
strategy (clr+pseudo) to around 0.9 for cumulative sum scaling normalization
with pseudo-zero replacement strategy (CSS+pseudo). This implies that for a
given taxonomic level, the resulted weighted F1 score will be highly
influenced by the normalization and zero replacement strategy. Traditionally,
microbiome researchers use the total sum scaling normalization without any
zero replacement strategy (TSS+none) on their data which has a range of 0.80
to 0.90 weighted F1 scores for the RF model (0.8 to 0.85 for the BNN model)
depending on the taxonomic level.

The strong interaction effects of taxonomic level, normalization and zero
replacement strategy prevent us from making recommendations about the best
data preprocessing practices that can be generalizable to other datasets. We
conclude by suggesting data practitioners to try different normalization and
zero replacement strategies rather than just one, but see Section B.2.5 for
more recommendations.

Figure 3: Weighted F1 scores (y-axis) for random forest and Bayesian neural
network (Bayesian NN) models for the pitted scab disease (Scabpit) under the
20 normalization/zero replacement strategies (x-axis). The lack of pattern
prevents us from making recommendations of optimal strategies for microbiome
data. We can conclude, however, that taxonomic levels, normalization and zero
replacement strategies have an effect on the prediction accuracy of the models
as evidenced by the broad range displayed by the points.

####  B.2.3 Effective preservation of predictive signal with different feature
selection strategies

One of the standard steps in the machine learning pipeline is feature
selection, especially for cases of high-dimensional data. We compare the
ability to retain predictive signal of three feature selection strategies:
standard importance score from machine learning methods, comparison of
microbial network topologies and combination of both. More details on the
feature selection strategies can be found in Methods.

Figure 4 shows the weighted F1 scores for the two types of models (RF and BNN)
on pitted scab disease (Scapbit) under different subsets of predictors: (1)
all OTUs (ALL-OTU), (2) only OTUs that were identified as important by the ML
strategy (OTU-S1), (3) only OTUs that were identified as important by the
network comparison strategy (OTU-S2), (4) OTUs that were identified as
important by both strategies (OTU-S3), or (5) OTUs that were not identified as
important by neither strategy (OTU-S0). For fair comparison, we include the
same number of predictors in OTU-S0 as in OTU-S3. Similar figures for other
responses are shown in Figures S11 to S15 in the Appendix.

Again, we perceive strong interaction between taxonomic level and feature
selection strategy. For the RF model, the highest weighted F1 score is
achieved when including all OTUs (ALL-OTU) at Order level whereas for the BNN
model, the highest weighted F1 score is achieved when including OTUs
identified by the ML strategy (OTU-S1) at Genus level. RF models on all OTUs
(ALL-OTU) have weighted F1 score above 0.8 in all taxonomic levels which
suggest that this model could be a better alternative compared to BNN which is
more computationally intensive. There are also smaller differences in RF
models when comparing the performance on OTU-S3 (important OTUs) and ALL-OTU
(all OTUs) which suggests that the feature selection strategy is sufficient to
preserve the predictive signal in the data while reducing the number of
predictors in the model. This is relevant for computationally intensive models
such as BNN that do not allow the inclusion of all OTUs for certain taxonomic
levels.

Figure 4: Weighted F1 scores (y-axis) by random forest and Bayesian neural
network (Bayesian NN) models for the pitted scab disease (Scabpit) by feature
selection strategy (x-axis) including all OTUs (All-OTU), OTUs selected by the
ML method (OTU-S1), the network comparison method (OTU-S2), both methods
(OTU-S3), or neither method (OTU-S0).

####  B.2.4 Robust prediction in imbalanced datasets with data augmentation

High prediction power in imbalanced datasets is misleading as a naive
predictor that classifies all samples as the majority class will have high
accuracy. In our data, black scurf disease is highly imbalanced, and thus, the
high prediction accuracy is unreliable. We confirmed, however, that after data
augmentation which balanced the data, accurate prediction persisted.

To illustrate this, Figure 5 depicts the weighted F1 scores on original and
augmented datasets for all yield and disease outcomes (x-axis) and both models
(RF and BNN). The range of each box plot depicts the weighted F1 scores for 20
normalized datasets at each taxonomic level. We observe that black scurf and
pitted scab can be reliably predicted across taxonomic levels as their median
weighted F1 scores for all taxonomic orders is around $0.8$ when the models
are fitted on the original datasets. As mentioned before, however, black scurf
is highly imbalanced, so the results on the original data are not reliable.
Fortunately, the median weighted F1 scores on augmented data (which is
perfectly balanced by design) increases for both diseases, such that they are
around $0.9$ for all taxonomic levels. These results suggest that data
augmentation, especially in cases of highly imbalanced data, is an appropriate
strategy that increases the accuracy of the prediction. One has to be careful,
however, in that augmented data can yield certain models prohibited. For
example, the BNN model could not be fit on the augmented datasets for Order,
Family nor Genus levels due to computational limitations.

Figure 5: Boxplots with weighted F1 scores by model (rows): random forest (RF)
and Bayesian neural network (BNN) for original and augmented OTU predictors
(columns). The x-axis represents the yield and disease outcomes. The range of
each box plot depicts the weighted F1 scores for 20 normalized datasets at
each taxonomic level.

####  B.2.5 Identifying general practices to predict potato disease from
microbiome data using full model selection models

As evidenced by our analyses, every single data and model choice has an impact
on the predictive performance of our methods. The effects of different data
preprocessing steps appear to strongly interact, and thus, we could not
identify clear patterns on strategies to maximize prediction power.

With a full model selection (FMS) [25, 26, 24] strategy, however, we are able
to identify the choices that yield the highest measure of performance. More
details on the FMS models can be found in Methods. Figures 6 and 7 show the
FMS decision trees for the RF and BNN models on pitted scab disease,
respectively. A FMS decision tree shows the different data preprocessing steps
that yield different weighted F1 scores, so that practitioners can select the
options that result in the highest predictive power. Here, we have five
taxonomic levels, 20 normalization+zero replacement strategies, and 2 data
augmentation options: no data augmentation (Aug=0) and data augmentation
(Aug=1). Thus, in total, we have 200 data preprocessing options (20
normalization strategies times 5 taxonomic levels times 2 data augmentation).

To interpret a FMS decision tree, each node corresponds to a specific step in
the data preprocessing pipeline, for example, whether to perform data
augmentation or not. If the condition is true, we follow the branch to the
left; if the condition is false, we follow the branch to the right. At the top
of the decision tree, we have the root which represents the data preprocessing
step that has the greatest effect on model accuracy. At the bottom of the
decision tree, we have the leaves with the average weighted F1 score of the
model fitted on the data that satisfy all conditions towards the root. The
darker the red, the higher the weighted F1 score. Each node also displays the
percentage of data preprocessing options included in the node. For example, in
Figure 6, the root node covers 100% of the options with average weighted F1
scores 0.865. The condition at the root node ($\text{Aug}=0$) represents the
case of "no data augmentation". Thus, "true" (left of the root) means "no data
augmentation", and "false" (right of the root) means "data augmentation". For
simplicity, we denote the 20 normalization/zero replacement strategies as
$\text{NM}_{i}$ for $i=1,\dots,20$. See Section D.1.2 for a description on
each normalization/zero replacement strategy.

For the FMS decision tree for the RF model (Figure 6), the highest weighted F1
score (0.934 with 0.5% of the data) is achieved with data augmentation,
normalization/zero replacement strategy #6 (CSS+pseudo), and Order level.
Another path of the decision tree follows data augmentation and any
normalization/zero replacement strategy except #6 (CSS+pseudo), #14
(rarefy+pseudo), and #18 (clr+pseudo) which yields an average weighted F1
scores of 0.892 for 42.5% of the data preprocessing options. If data
augmentation is not an option (left of the root), the highest weighted F1
score available is 0.868 with Phylum level, and any normalization/zero
replacement strategy except #18 (clr+pseudo) or #20 (clr+bayesMult). For the
FMS decision trees on the other responses, see Figures S16 to S20 in the
Appendix.

Similarly, in Figure 7 for the BNN model, the highest weighted F1 score (0.896
with 15% of the data preprocessing options) is achieved when we do data
augmentation, we use any taxonomic level except Phylum, and we use any
normalization/zero replacement strategy except #10 (COM+pseudo) and #18
(clr+pseudo). See Figures S21 to S25 in the Appendix for other responses.

While the specific recommendations on normalization, zero replacement and
taxonomic level are model-specific, both models perform better with data
augmentation. In terms of taxonomic level, we note that the BNN was only run
on Phylum and Class, and thus, the highest accuracy is obtained with Class
level (when Phylum$=0$ is true). This does not contradict the result from the
RF that identified Order level as the one yielding higher accuracy. We cannot
rule out that the BNN would also have higher accuracy with Order compared to
Class. The results from the RF, though, seem to suggest that there is a peak
at Order, and more granularity in Family and Genus does not seem to provide
more predictive power.

Figure 6: Full model selection decision tree with a maximum depth of 4
summarizing the results of random forest models on pitted scab disease
(Scabpit). When the condition at a node is true, we follow the branch on the
left, and when the condition is false, we follow the branch on the right. The
percentage of the data preprocessing options and mean of weighted F1 scores
are shown in each nodes. The darker the red, the higher the weighted F1 score.
Figure 7: Full model selection decision tree with a maximum depth of 4
summarized the results of Bayesian neural network models on pitted scab
disease (Scabpit) disease. When the condition at a node is true, we follow the
branch on the left, and when the condition is false, we follow the branch on
the right. The percentage of the data preprocessing options and mean of
weighted F1 scores are shown in each nodes. The darker the red, the higher the
weighted F1 score.

###  B.3 Prediction of Potato Disease from Environmental Data

One of the questions to address in our work is whether prediction accuracy is
improved by the inclusion of microbiome data, or if environmental factors
provide enough signal to classify potatoes in diseased or non-diseased groups.

We found that environmental factors indeed provide sufficient signal to
predict pitted scab disease as illustrated in Figure 8 which shows the
weighted F1 scores by RF and BNN models based on environmental (soil
characteristics) data for pitted scab. The range of each boxplot corresponds
to the six scaling methods described in Section D.1.2. In contrast with the
normalization methods in microbiome data, we observe here that the scaling
methods do not seem to have an effect on prediction as evidenced by narrow
boxplots, and that weighted F1 scores are all higher than 0.7, and therefore,
comparable to the models fitted on microbiome data alone. These results
suggest that environmental factors alone are powerful to predict the incidence
of pitted scab in the tubers. As microbiome data is more expensive than
environmental data, we suggest to prefer environmental predictors under
restricted monetary budget. See Figures S26 to S30 in the Appendix for other
responses.

Figure 8: Boxplots of the weighted F1 scores by random forest (RF) and
Bayesian neural network (BNN) models for pitted scab disease (Scabpit) using
environmental predictors.

###  B.4 Leveraging Microbiome and Environmental Data in the Prediction of
Potato Disease

As expected, prediction accuracy is improved by the inclusion of both
datatypes: microbiome and environmental. Figure 9 shows the weighted F1 scores
by RF and BNN models based on combined datasets with environmental and
microbial predictors for pitted scab. We only focus on the most accurate
models identified in Section B.1. First, we note that a model that uses OTU
abundances outperforms a model that uses alpha diversity as a predictor
(comparison of Alpha with OTU-S3) for both types of models (RF and BNN). This
suggests that we lose information by transforming abundances into diversity
measures. Second, models including only OTU abundances (OTU-S3) perform
comparably to models that include both types of predictors (OTU-S3+Soil+DS)
which suggests that the microbial data indeed has substantial predictive power
on its own, but adding microbiome to soil predictors may not provide much
benefit for high predictive power, with the only exception of Phylum level
OTU-S3+Soil+DS in a RF model (Figure 9). Generally, the model with only soil
information (shown as a blue dashed line) performs just as accurately. Third,
contrary to prior expectations that microbial communities at finer resolution
would be a better choice for predicting pitted scab or other diseases, our
study does not find any evidence that the prediction power increases when
moving up from Phylum to Genus level. Particularly, the prediction power of
OTU-S3 in RF model increases from Class to Genus, and this pattern is not
preserved when diversity is used instead of OTU abundances. For example, a
model with only alpha diversity as the predictor (Alpha) shows decreasing
weighted F1 score as we move from Phylum to Genus level. Both models (RF and
BNN) when including all types of predictors (OTU-S3+Soil+DS) result in the
similar weighted F1 score regardless of taxonomic level. See Figures S31 to
S35 in the Appendix for other responses.

Figure 9: Boxplots with the weighted F1 scores (y-axis) by random forest and
Bayesian neural network (Bayesian NN) models (x-axis) for pitted scab disease
(Scabpit). The range of each box plot depicts the weighted F1 scores for
normalized datasets at each taxonomic level. There are 6 normalization methods
for alpha diversity and environmental factors (Soil and DS), and one
normalization method for microbiome (OTU-S3). See Table 3 and 4 for
description of models and number of normalization methods for each predictors.
The models including both types of predictors outperform other models, yet
models including microbiome data alone (OTU-S3) are comparable which suggests
that the microbial information indeed contains signal to predict the disease
outcome on its own. However, microbiome data is more expensive to collect, and
perhaps not necessary, given that the model without microbiome data (Soil)
performs just as accurately (blue dashed line).

##  Appendix C Discussion

Soil microbiome represents the most complex and least understood aspect of
soil health. In this study, we use ML techniques such as random forest (RF)
and Bayesian neural network (BNN) to determine whether soil microbial
information has any predictive power for plant outcomes such as yield and
disease. Our results show that the microbiome plays a role in predicting
potato disease, but the most accurate prediction originates from a combination
of both datatypes: microbiome and environmental factors. Given that prediction
with environmental factors alone was sufficiently powerful, it is uncertain
whether the extra expense to sequence microbiome data is worth the cost.

##### Can machines classify what humans cannot? The importance of accurate
labels.

The prediction power of microbiome data varies depending the outcome we want
to predict. For example, among all models that predict diseases, models for
the pitted scab disease receive very high weighted F1 scores compared to other
diseases. We further confirm this predictive power by comparing the
performance of models trained on the real microbiome data to models trained on
randomly generated data. Given that the prediction of the pitted scab disease
is far from random, we can confidently conclude that this disease can be
accurately predicted from microbiome data. It is noteworthy, however, that
pitted scab disease is precisely one of the diseases that are easier to be
visually detected, and thus, there is a reliable separation among the two
classes (diseased and non-diseased) which is aiding in prediction by machine-
learning models. Other diseases, and more so yield, do not have such clear
distinction between classes which results in lower predictive power. That is,
we believe that the lack of prediction accuracy in yield is not driven by a
lack of a biological connection between soil microbiome and yield, but on the
lack of accurate labels that distinguish the two classes (e.g. low and high
yield). If humans cannot distinguish what is low yield vs high yield, then
that ambivalence will propagate into the machine-learning classification. This
conclusion seems to be confirmed when we notice that yield cannot be
accurately predicted by any of 14 models in consideration (Table 3). We
conclude that one of the main challenges when applying ML methods in
biological applications is the artificial binarization of phenotypes. More
work is needed to improve the performance of regression models that can
predict continuous phenotypes when faced with limited sample sizes that are
common in biological domains.

##### Data preprocessing has a substantial impact in prediction performance.

We show that the normalization and zero replacement strategies for microbiome
data sets have a huge impact on results, and unfortunately, we could not find
any pattern to suggest that one strategy is superior to others. The full model
selection (FMS) decision tree, however, suggests that certain normalization
methods such as centered log-ratio normalization with pseudo-zero replacement
(clr+pseudo) or centered log-ratio normalization with Bayesian-multiplicative
treatment (clr+bayesMult) should not be used for disease prediction for future
studies. We recommend domain scientists to perform different types of
normalization and zero replacement strategies in their data prior to be used
as input for prediction. See also [4, 22, 7] for other performance comparisons
on types of normalizations.

In addition, the data augmentation proved beneficial for our models in various
ways. First, it added credibility for imbalanced responses such as the black
scurf disease, as adding augmented data to the training set enables the model
to train on both response labels of equal size, preventing the model from
making biased decisions towards one label. Second, the model offers slightly
better performance when training on the augmented data for all disease
responses. Finally, data augmentation yields a more robust model, which lowers
the chance of overfitting on the current data and should offer better
performance on future data if we want to make decisions based on our current
models. We advice practitioners, however, to be careful with data augmentation
techniques to prevent data leakage. Unlike data augmentation in imaging data
which is straight-forward, data augmentation of relative abundances can
artificially boost accuracy by having the original sample in the training set
and the augmented sample in the testing set. More work is needed to produce
reliable data augmentation techniques for biological data that are not images
beyond the standard Gaussian noise incorporation.

##### Feature selection effectively preserves the predictive signal on lower
dimensions.

In terms of feature selection, we considered different strategies to select
important features (ML, network comparison, and intersection of both). There
is significant overlap between important OTUs by two methods and the inclusion
of this subset of predictors allowed us to build less complex models with
comparable good performance (Figure 4 for pitted scab disease). We did not
observe one feature selection strategy that outperformed the others. For
example, the weighted F1 scores obtained when using the OTUs identified as
important by the ML methods are comparable to the weighted F1 scores when
including OTUs identified as important by network comparison. It is worth
mentioning, however, that the performance is also comparable to that obtained
when including all OTUs which shows that the feature selection strategies work
at preserving the OTUs that have predictive signal while simultaneously
allowing computationally expensive models (like BNN) to be applied. Tables S2
and S3 in the Appendix contain the important OTUs identified by the two
feature selection methods for the Phylum level associated with pitted scab
disease.

##### Finer taxonomic levels do not provide higher predicting power in all
scenarios.

In our study, we could not identify an overall prediction superiority on
certain taxonomic levels compared to others. While it is intuitive to expect
that finer taxonomic levels would provide more power to predict potato
outcomes, this is only true in some cases such as the OTU-S3 model in RF and
the OTU-S1 model in BNN (Figure 4) where Genus level indeed outperforms the
other taxonomic levels. However, we cannot find evidence to support this
expectation in all other cases.

Even though this finding disagrees with common expectations in microbiome
research, we believe that it is a positive result from the machine learning
perspective. If less refined taxonomic levels perform just as well as finer
taxonomic levels, then the model could be fitted with fewer features, and
would thus be less computational expensive and more efficient. Furthermore,
the more features in the data, the more samples are needed to train the
models. So, ML models that are trained on limited samples sizes (as the ones
we have in this study) would benefit greatly by having fewer features.

##### Limited predictive power in soil microbiome compared to environment.

When including environmental features such as soil physicochemical properties
and microbial population density of soil in the model, we achieved higher
weighted F1 score values. In general, we investigated 14 different models for
yield and disease prediction with different combinations of microbiome and
environmental data. Results show poor performance in predicting yield across
different models. The best results for pitted scab are achieved by combining
alpha diversity and soil chemistry (Alpha+Soil) for RF and important OTUs and
soil data (OTU-S3+Soil) for BNN. Although the best-performing models include
microbiome predictors (Alpha and OTU-S3), it is important to note that the
models without microbiome data are comparably powerful as those including
microbiome data. While the models trained with microbiome data alone show that
microbiome can effectively predict pitted scab disease, the fact that models
without this type of data continue to perform well provides evidence that
microbiome data may not be necessary to achieve reasonable prediction. In
fact, when collection of microbiome data requires much higher cost investment,
the increased prediction accuracy by including microbiome predictors may not
be enough to justify the extra cost.

The RF and BNN methods offer similar performance in all models, which
validates the results of both methods. However, in this study, BNN take far
longer to train due to weight sampling and approximation, especially on more
refined taxonomic level such as Order, Family, and Genus. Thus, we believe
that based on the current dataset, RF is the preferred decision-making model
with high prediction potential for disease outcomes when including a
combination of microbiome and environmental predictors.

##  Appendix D Methods

Figure 1 shows a graphical representation of our pipeline with three major
steps: (1) data preparation, (2) feature selection, and (3) classification
based on random forest (RF) and Bayesian neural network (BNN) models with
various types of predictors, including microbiome data (OTUs), environmental,
and a combination of both. We describe each step in the pipeline in the next
subsections.

###  D.1 Data description

In this study, we focus on the soil microbiome (matrix of abundances) in a
variety of taxonomic orders, including Phylum, Class, Order, Family, and Genus
as well as other environmental information from soil samples acquired from
potato fields in Wisconsin and Minnesota. The dataset consists of measurements
related to soil health, potato yield and soil quality information. The soil
health data are collected from pre-plant fields and include soil
physicochemical properties, soil microbiome composition, soil microbiome
diversity, and soil pathogen abundance. The potato yield and quality
assessments are taken from the following growing season with indicators
including tuber yield and disease severity. Overall, we have 256 samples, 108
of which are taken from fields in Minnesota, and 148 of which are taken from
fields in Wisconsin. We list all measurements in Table S1 in the Appendix.

####  D.1.1 Binarization of Response Variables

We have six phenotypes (response variables) of interest, four of them
correspond to diseases and two of them correspond to yield (Table 1).

Response | Description
---|---
Yield_Meter | Weight of potatoes (grams) per one meter yielded at harvest time.
Yield_Plant | Weight of potatoes (grams) per plant yielded at harvest time
Scab | Number of tubers that carry scab (superficial+pitted) disease on the sampled plant.
Scabpit | Number of tubers that carry pitted scab disease on the sampled plant.
Scabsuper | Number of tubers that carry superficial scab disease on the sampled plant.
Black_Scurf | Percentage coverage of black scurf disease on the sampled plant
Table 1: Description of responses of interest in this study.

All six responses in the dataset are continuous, so we need to binarize them
to fit the classification models. For the disease-related responses, we simply
make the binary label $0$ if there is no presence of disease, and $1$ if there
was detection of disease (that is, if the continuous response is greater than
$0.0$). Binarizing the yield responses variables is harder as there is not a
universal standard to classify potato yield to be low or high. Furthermore,
yield values are highly dependent on the type of potato variety. We assign a
label of $0$ (low yield) to samples with a yield less than the variety-
specific median. Similarly, we assign a label of $1$ (high yield) to samples
with yield greater than the variety-specific median. We illustrate this
approach in Figure S1 in the Appendix.

After binarization, we note that pitted scab disease (denoted Scabpit in the
figures), superficial scab disease (denoted Scabsuper in the figures), and
both yield responses are balanced, whereas other responses are highly
imbalanced: scab disease (denoted Scab in the figures) has 80% of samples
labeled as $1$, and black scurf disease (denoted Black_Scurf in the figures)
has only six samples labeled as $1$. We use these imbalanced cases to assess
the performance of the methods under imbalance settings and data augmentation
strategies.

####  D.1.2 Data Filtering, Normalization, and Zero Replacement

The input data is a matrix with non-negative read counts that were generated
by a sequencing procedure. Let $w^{(k)}=[w_{1}^{(k)},...,w_{p}^{(k)}]$ be the
total read counts of sample $k$ containing $p$ OTUs, where $w^{(k)}$ is a
composition that adds up to a fixed value of
$m^{(k)}=\sum_{i=1}^{p}w_{i}^{(k)}$. This value $m^{(k)}$ is the sequencing
depth, which varies across samples and is predetermined by technical factors
resulting in highly sparse data. It is reasonable to filter out a certain set
of OTUs as first data preparation step. For filtering, we only include OTUs
that appear in at least 15 samples.

Table 2 displays the number of features (OTUs) before and after filtering for
different taxonomic levels.

Level | # OTUs | # OTUs after filtering
---|---|---
Phylum | 57 | 42
Class | 152 | 108
Order | 378 | 224
Family | 476 | 255
Genus | 1319 | 485
Table 2: Number of OTUs per taxonomic level (first column) in the original
data (second column) and after filtering out OTUs that do not appear in at
least 15 samples (third column).

As mentioned, the input data is compositional and highly sparse. It is known
that ML methods do not perform well with unnormalized data [18] and with
sparse data [29]. Therefore, we explore the effect of four zero replacement
strategies (to overcome sparsity) and five normalization strategies (to
overcome compositionality). All strategies are implemented in the NetCoMi R
package [23].

In particular, we consider the four zero replacement strategies: (1) the
original dataset which included zeros (denoted none in the figures), (2)
pseudo-zero replacement which replaces zero counts by a predefined pseudo
count (denoted pseudo in the figures), (3) multiplicative zero replacement
which imputes left-censored compositional values by a given fraction and
applies a multiplicative adjustment to preserve the multivariate compositional
properties of the samples (denoted multRepl in the figures) [16], and (4)
Bayesian-multiplicative treatment which imputes zero counts by posterior
estimates of the multinomial probabilities generating the counts, assuming a
Dirichlet prior distribution (denoted bayesMult in the figures) [17].

Next, we use five normalization methods: (1) Total sum scaling which simply
converts counts to appropriately scaled ratios (denoted TSS in the figures)
[4], (2) Cumulative sum scaling which rescales the samples based on a subset
(quartile) of lower abundant taxa, thereby excluding the impact of highly
abundant taxa (denoted CSS in the figures) [4], (3) Common sum scaling in
which counts are scaled to the minimum depth of each sample (denoted COM in
the figures) [18], (4) Rarefying which random samples without replacement
after a minimum count threshold has been applied (denoted rarefy in the
figures) [10], and (5) Centered Log-ratio which transforms the data using the
geometric mean as the reference (denoted clr in the figures) [2].

With four zero replacement methods and five normalization methods, we create
20 datasets by the combination of zero replacement and normalization
techniques. Our goal is to study the effect of the zero replacement and
normalization choice in the performance of the deep learning methods. Namely,
we have the following 20 combinations, $\text{NM}_{1}$: TSS+none,
$\text{NM}_{2}$: TSS+pseudo, $\text{NM}_{3}$: TSS+multRepl, $\text{NM}_{4}$:
TSS+bayesMult, $\text{NM}_{5}$: CSS+none, $\text{NM}_{6}$: CSS+pseudo,
$\text{NM}_{7}$: CSS+multRepl, $\text{NM}_{8}$: CSS+bayesMult,
$\text{NM}_{9}$: COM+none, $\text{NM}_{10}$: COM+pseudo, $\text{NM}_{11}$:
COM+multRepl, $\text{NM}_{12}$: COM+bayesMult, $\text{NM}_{13}$: rarefy+none,
$\text{NM}_{14}$: rarefy+pseudo, $\text{NM}_{15}$: rarefy+multRepl,
$\text{NM}_{16}$: rarefy+bayesMult, $\text{NM}_{17}$: clr+none,
$\text{NM}_{18}$: clr+pseudo,, $\text{NM}_{19}$: clr+multRepl, and
$\text{NM}_{20}$: clr+bayesMult. For convenience, we use the notation
$\text{NM}_{i}$ (Normalization Method) for $i=1,\dots,20$ in the Full Model
Selection section (Section B.2.5).

For the environmental predictors of soil chemistry and microbial population
density in the soil, we apply six scaling methods: (1) standardize features by
subtracting the mean and scaling to unit variance [8]; (2) scale each feature
to a $[0,1]$ range; (3) scale each feature by its maximum absolute value; (4)
scale features by subtracting the median and scaling to the interquantile
range [32]; (5) transform the features to follow a normal distribution [21];
(6) normalize samples individually to the unit norm.

After normalization, the datasets are split into training, validation, and
testing sets with 10-fold cross-validation. We used 80% of samples for
training and validation, and 20% for testing.

####  D.1.3 Data Augmentation

There are three main goals that we wish to achieve with data augmentation: (1)
improve the model’s prediction performance with more artificial samples; (2)
balance the number of labels with artificial samples, and (3) make the model
more robust and avoid overfitting with unseen (artificial) data. We note that
augmenting the whole dataset and then splitting into training and testing sets
would result in data leakage. For example, when the original sample is in the
testing set and the augmented sample from this sample is in the training set,
the model is essentially training and testing on the same sample since the
normalized values of OTUs are very close. Thus, we split the data into
training and testing sets first and only augment the training set. This
strategy also allows us to have a fair performance comparison for augmented
and non-augmented sets with the same testing data.

Regarding the data augmentation procedure, instead of simply adding a randomly
generated noise to the original sample, we subset the data by variety and
label, compute the mean (and standard deviation) abundance value for this
subset, and create a new sample that includes the original data plus a
Gaussian error with mean $\mu/100$ and standard deviation $\sigma/100$ where
$\mu,\sigma$ are the subset-specific mean and standard deviation,
respectively. This approach is illustrated in Figure S2 in the Appendix. By
the end of this procedure we would have a balanced augmented training set with
400 samples per label for each of the five taxonomic levels (Phylum, Class,
Order, Family, and Genus).

###  D.2 Feature Selection

Feature selection involves the identification of important features (or
covariates) that have high predictive power. Given the high-dimensionality of
the data (e.g. 256 original samples for 485 OTUs at the Genus level), feature
selection is necessary, especially for BNN models that are computationally
expensive.

We pursue two approaches for feature selection: (1) using machine-learning
models to assess variable importance, and (2) using network analyses. To focus
exclusively on the effect of feature selection, we only consider one type of
normalization and zero replacement strategy in this investigation, namely,
total sum scaling normalization without zero replacement ($\text{NM}_{1}$:
TSS+none).

####  D.2.1  Using ML Models for Feature Selection

To identify important OTUs, we use six ML strategies implemented in scikit-
learn [21]: (1) "SelectKBest" method selects features based on the $k$ highest
analysis of variance F-value scores, (2) select the top $k$ features based on
the mutual information statistic, (3) recursive feature elimination (RFE) with
logistic regression, (4) RFE with decision tree, (5) RFE with gradient
boosting, and (6) RFE with RF. In addition to the six ML strategies, we
consider a 7th strategy which consists in including OTUs in the model if their
maximum value is within the top 30%.

After running all seven strategies, we assign a value ("TOTAL") to each OTU
based on the number of times the OTU is selected as an important feature under
the seven criteria. That is, an OTU that is selected as important by all seven
strategies will have a value of 7. The OTUs are sorted based on "TOTAL" column
and top 30% of them are selected as important features. Thus, 30, 36, 75, 85,
and 162 OTUs are selected for Phylum, Class, Order, Family and Genus levels,
respectively. For example, Table S2 in the Appendix shows the top 30 OTUs at
Phylum level and by which strategies they are identified as important features
to predict the pitted scab response.

####  D.2.2 Using Network Comparison for Feature Selection

Next, we identify important OTUs by comparing their interactions in microbial
networks when the network is constructed with samples from one class (say, low
yield) versus with samples from the other class (say, high yield). Indeed,
evidence shows that microbial interactions can help differentiate between crop
disease states [23]. We identify the OTUs that interact differently on the
samples corresponding to one class versus another class as relevant OTUs that
contain information about the outcome.

We use the SPRING method (semi-parametric rank-based approach for inference in
graphical model) [33] in the NetComi package to construct two microbial
networks: one corresponding to label $0$ and one corresponding to label $1$.
To build these networks, the estimated partial correlations are transformed
into dissimilarities via the signed distance metric and the corresponding
similarities are then used as edge weights. We then compare these networks in
order to identify OTUs that have the greatest difference based on degree
values. Figure 10 shows an example of two microbial networks with different
graphical structures, and Table S3 in the Appendix lists OTUs ranked by the
difference in degree in the two microbial networks for diseased and non-
diseased classes of pitted scab response. The OTUs are sorted based on the
values of degree difference column and similar to ML strategy top 30% of them
are selected as important features.

Figure 10: Microbial networks at Phylum level for non-diseased (Label $0$) and
diseased (Label $1$) samples for pitted scab (Scabpit). We use the SPRING
method [33] to reconstruct these microbial networks. Node colors represent
clusters, which are determined using greedy modularity optimization. Clusters
have the same color in both networks if they share at least two OTUs. Green
edges correspond to positive associations and red edges to negative ones.
Nodes that are unconnected in both groups are removed. OTUs names are
abbreviated (see Table S3 in the Appendix for the original names).

####  D.2.3 Combination of ML and Network Strategies in Feature Selection

We define a scoring value for each OTU based on whether they are identified as
important by ML strategy ($score=1$), by the network comparison ($score=2$),
or both ($score=3$). If the OTU is not identified as important by any
strategy, it is denoted $score=0$. Figure 11 displays the scoring values of
all OTUs at the Phylum level for all responses.

Figure 11: Heatmap with scoring of OTUs (Phylum level) as important for each
response. Y-axis shows the responses and x-axis shows the OTUs in Phylum
level. We define a scoring value for each OTU based on the selection by ML
strategy ($score=1$), network comparison strategy ($score=2$), both
($score=3$) or neither ($score=0$).

###  D.3 Model descriptions: random forest and Bayesian neural network

We apply two types of classification models to predict potato disease and
yield: random forest (RF) and Bayesian neural networks (BNN). We consider
separately three types of predictors: OTU abundances (5 taxonomic levels and
20 normalization/zero replacement strategies described in Section D.1.2),
environmental predictors such as soil characteristics and microbial population
density, and a combination of both types. Table 3 lists all the models we
consider in this study.

Type of predictors | Name of the model | Predictors included in the model
---|---|---
OTU | ALL-OTU | OTU abundances
OTU-S0 | OTUs with a score of zero (not selected by ML or network comparison
|  | feature selection strategies)
| OTU-S1 | OTUs selected by the ML feature selection strategy (score of one)
| OTU-S2 | OTUs selected by the network comparison feature selection strategy
|  | (score of two)
| OTU-S3 | OTUs selected by the ML and network comparison feature selection
|  | strategies (score of three)
| Alpha | Alpha diversity
Environmental | Soil | Soil chemistry
DS | Microbial population density in soil
| Soil+DS | Combination of soil chemistry and microbial population density in soil
Combination | Alpha+Soil | Combination of alpha diversity and soil chemistry
Alpha+Soil+DS | Combination of alpha diversity, soil chemistry, and
|  | microbial population density in soil
| OTU-S3+Soil | Combination of OTUs with score of three and soil chemistry
| OTU-S3+DS | Combination of OTUs with score of three and microbial population
|  | density in soil
| OTU-S3+Soil+DS | Combination of OTUs with score of three, soil chemistry, and
|  | microbial population density in soil
Table 3: Models under study for three types of predictors (first column):
OTUs, environmental and both. The name of each model (second column) is used
in figures and tables in the text.

Table 4 shows the number of predictors included in each model as well as the
data choices related to taxonomic level or normalization and zero replacement
strategies. For example, using all OTUs (first row), we have five taxonomic
levels and 20 normalization+zero replacement strategies each, so in total, we
have 100 (20 normalization strategies times 5 taxonomic levels) different OTU
datasets. For each of these 100 datasets, the number of predictors (i.e.,
OTUs) would depend on the taxonomic level being analyzed. For example, at the
Phylum level, there might be 42 predictors, while at the Genus level, there
could be up to 485 predictors.

Model | Taxonomic levels | Normalization/zero replacement strategies | Number of predictors
---|---|---|---
ALL-OTU | 5 | 20 | 42–485
OTU-S0 | 5 | 1 | 7–54
OTU-S1 | 5 | 1 | 5–134
OTU-S2 | 5 | 1 | 5–134
OTU-S3 | 5 | 1 | 7–54
Alpha | 5 | 6 | 9
Soil | – | 6 | 12
DS | – | 6 | 4
Table 4: Types of predictors per model, number of taxonomic levels,
normalization/zero replacement strategies, and range/number of predictors for
each type.

####  D.3.1 Random Forest Model

The RF classifier is a powerful ML technique that has gained significant
popularity in the last two decades because of its accuracy and speed. RF
randomly creates an ensemble of decision trees. Each tree picks a random set
of samples (bagging) from the data and models the samples independently from
other trees. Instead of relying on a single learning model, RF builds a
collection of decision models, and the final decision is based on the output
of all the trees in the model. The bagging approach promotes the generation of
uncorrelated trees which reduces the risk of overfitting.

Each decision tree is generated individually without any pruning and each node
is split using a user-defined number of features. By expanding the forest to a
user-specified size, the technique generates trees with a high variance and
low bias. The final classification choice is determined by summing the class-
assignment probability obtained by each tree. A new unlabeled data in testing
set input is thus compared to all decision trees formed in the ensemble, with
each tree voting for class membership and the membership category with the
most votes will be picked.

The RF has several hyper parameters to be determined by the user, such as the
number of decision trees to be generated, the number of variables to be
selected and tested for the best split when growing the trees, the maximum
depth of the tree, the minimum number of samples required to split an internal
node, among others. Generally, a grid search is combined with K-fold cross
validation to select the best hyper parameters [30]. GridsearchCV is a well-
known search method which is available in scikit-learn [21] and it evaluates
all possible parameter combinations to determine optimal values.

Here, we set different values for parameters (see Table S4 in the Appendix)
and tune them using GridsearchCV to find the optimal values for the RF
classifier. GridSearchCV uses a “score” method for evaluating the performance
of the cross-validated model on the test set. We use the F1 score to evaluate
the performance of the model. The F1 score is the harmonic mean of precision
and recall, where the proportional contributions of precision and recall to
the F1 score are equal. A F1 score of 1 is the best, while a score of 0 is the
worst. Finally, when the parameters for the RF model are tuned, we use 20% of
samples to report the performance of the final model. We use the weighted
average F1 score to evaluate the performance of the models, which is computed
by averaging all the per-class F1 scores while accounting for the number of
samples in each class.

####  D.3.2 Bayesian Neural Network Model

Given that neural network models are intrinsically data-hungry and microbiome
applications have relatively small sample sizes, BNN models provide a suitable
alternative for small sample sizes as they provide natural protection against
overfitting. This is due to the fact that distributions are considered for the
parameters in the model which allow us to marginalize them so that the
prediction is based on data points alone [9].

We provide details on the BNN model formulation next. Let
$\mathcal{D}=\\{x_{i},y_{i}\\}_{i=1}^{n}$ be the data with $n$ samples where
$x_{i}$ is a $N$-dimensional feature vector in $\mathbb{R}^{N}$ and
$y_{i}\in\mathbb{R}$ is the target response. A BNN model has $L$ layers with
$K_{l}$ neurons in the $l^{th}$ layer. The weight parameters are defined in a
set of matrices $\mathcal{W}=\\{W_{l}\\}_{l=1}^{L}$ where each $W_{l}$
corresponds to the weights for the $l^{th}$ layer, and is of size
$(K_{l-1}+1)\times K_{l}$. For an input $x_{i}\in\mathbb{R}^{N}$, the neural
network maps it to a response $f(\mathcal{W},\textbf{x})$ by multiplying the
input of each layer by the weights and then transforming via an activation
function such as rectified linear unit (ReLU): $h(a)=\max(0,a)$.

Unlike regular neural network models, a BNN imposes a prior distribution on
the weights $\mathcal{W}\sim p(\mathcal{W})$ to capture their uncertainty. The
posterior distribution of the weights given the data is then given by:

$\displaystyle
p(\mathcal{W}|\mathcal{D})=\frac{p(\mathcal{W})\prod_{i=1}^{n}p(y_{i}|f(\mathcal{W},x_{i}))}{p(\mathcal{D})}$

which can be used to predict the response y* of unseen data x* via
$p(\textbf{y*}|\textbf{x*})=\int
p(\textbf{y*}|f(\mathcal{W},\textbf{x*}))p(\mathcal{W}|\mathcal{D})d\mathcal{W}$.

We use the original implementation of BNNs denoted “software for flexible
Bayesian modeling and Markov chain sampling" [20] for training and testing on
our data. For each BNN model, we use five hidden layers, one input layer and
one output layer. The input layer has $n$ input neurons with $n$ being the
number of predictors for the model. This number varies depending on the
predictors to include. For example, for the ALL-OTU model for the Phylum
level, we have 42 predictors (Table 2), and thus the input layer has $N=42$
neurons. For each hidden layer, we have $3N$ hidden neurons and one bias node.
For the output layer, we only have one neuron with a binary value for the
response. We use the hyperbolic tangent as the activation function:
$f(x)=tanh(x)$.

We choose a hierarchical zero-mean Gaussian prior on all weights and on bias
nodes. The prior for the weight from node $i$ of one layer to node $j$ of
another layer is formulated as follows [20]:

$\displaystyle\mathcal{P}(w_{ij}|\sigma_{w_{i}},\sigma_{a_{j}})=\frac{1}{\sqrt{2\pi}\sigma_{w_{i}}\sigma_{a_{j}}}\exp(-\frac{w_{ij}^{2}}{2\sigma_{w_{i}}^{2}\sigma_{a_{j}}^{2}})$

where $\sigma_{w_{i}}$ are the hyperparameters for the standard deviation of
the priors for node $i$ of layer $w$, and $\sigma_{a_{j}}$ is the hyper
parameter for the standard deviation of the priors for node $j$ of layer $a$.
Let $\tau_{w_{i}}=\sigma_{w_{i}}^{-2}$, then the hyper prior is given by

$\mathcal{P}(\tau_{w_{i}}|\tau_{w})=\frac{(\alpha_{w_{i}}/2\tau_{w})^{\alpha_{w_{i}}/2-1}}{\Gamma(\alpha_{w_{i}}/2)}\tau_{w_{i}}^{\alpha_{w_{i}}/2-1}\exp(-\frac{\tau_{w_{i}}\alpha_{w_{i}}}{2\tau_{w}})$
(1)

where $\alpha_{w_{i}}$ is the shape hyper parameter, and $\tau_{w}$ is the
scale hyper parameter. For weights between hidden layers and bias nodes, we
set $\alpha_{w_{i}}=2$ and $\tau_{w}=0.5$. For weights from the input layer to
the first hidden layer, we set a prior for $\tau_{w}$ as an inverse gamma
distribution (same as Equation 1) with shape hyper parameter $\alpha=2$ and
scale hyper parameter $\tau=0.5$. This extra hyper parameter diminishes the
effect of input nodes that have insufficient predictive power by setting
weights arbitrarily close to zero.

The training and the approximation of the posterior distribution are done via
a Hamiltonian (Hybrid) Monte Carlo (HMC) implemented in [20]. We select
different leap frog lengths and step sizes for different models. Leapfrog
lengths begin at 100 as starting point and are gradually decreased as a
leapfrog length of 100 is expected to deal optimally with challenging
estimation problems. Step size, on the other hand, is set at 0.1 as a starting
point and it is decreased and increased in search for an average rejection
rate smaller than 0.3 [19]. Burn-in is selected as the first half of the chain
so that weights are sampled on the second half of the converged chain.

Given that HMC does not scale well for high dimensional parameter spaces and
large datasets [31], we did not fit the BNN on all the models described in
Table 3. In particular, we do not fit a BNN model on Genus level for the ALL-
OTUs as this model would involve over 9 million weights in the network with
485 input neurons.

###  D.4 Predictive Power of Microbiome Data: Comparison to Random Data

One of the key questions we want to address in this study is whether the
microbiome data has predictive power to infer disease or yield outcomes. To
address this question more formally, we devise a simulation scheme in which we
simulate random microbiome data to train our RF model and compare the
prediction potential with the trained RF model on real microbiome data. If
there is indeed a signal to predict yield or disease outcomes within the real
data, then we expect the model trained with the real data to dramatically
outperform the model trained with random microbiome datasets. We use four
strategies to generate random datasets:

Strategy 1:

Random matrix with values between 0 and 1, normalized so that each row has a
sum of 1. We use the response vector from the original dataset.

Strategy 2:

Real microbiome matrix where each entry that is greater than zero is replaced
by a random number between 0 and 1. This strategy preserves the sparsity in
the real data as zero entries remain zero. Rows are normalized so that the sum
equals 1. We use the response vector from the original dataset.

Strategy 3:

Real microbiome matrix is used with a permuted response vector.

Strategy 4:

Rows in the real microbiome matrix are permuted, and the response vector is
unchanged.

For each strategy, we generate $N=200$ random datasets to train the RF model.
Let $F_{i}$ for $i=1,...,N$, denote the weighted F1 scores for the randomized
data and let $F_{original}$ be the weighted F1 score on the real data. Let
$X=\sum_{i=1}^{N}\mathbb{I}(F_{i}>F_{original})$ be the number of random
datasets that perform better than the read data. We use the exceeding value
(EV) $EV=\frac{X}{N}$ as test statistic to test the null hypothesis that the
real data performs just as random data in the prediction on disease outcomes.
The EV denotes the percentage of random datasets that perform better than the
real microbiome data out of $N$ random datasets. If $EV<\alpha$, we reject the
null hypothesis and conclude that there is predictive power (beyond random) on
the real microbiome data. Here, we set $\alpha=0.05$ as in [1].

###  D.5 Full Model Selection

Full model selection (FMS) [25, 26, 24] involves the process of listing all
data preprocessing steps, model options and selection of predictors, and using
a decision tree model to identify the choices that yield the highest measure
of performance. Here, we fit a FMS strategy with the following options: 1)
type of normalization, 2) type of zero replacement, 3) taxonomic level, and 4)
data augmentation. We combine the type of normalization and type of zero
replacement strategy into one variable (denoted $\text{NM}_{i}$ for
$i=1,\dots,20$). We focus on the weighted F1 score as measure of performance,
and we include all OTU predictors (that is, we do not consider feature
selection as one of the options to compare). We build the regression decision
tree by using the DecisionTreeRegressor which is available in scikit-learn
[21]. We use the default parameters in the DecisionTreeRegressor such as
"squared error" as the criterion to measure the quality of a split, a minimum
number of 2 for the samples required to split an internal node, and a minimum
number of 1 sample required to be at a leaf node. In order to create an
informative decision tree that can be interpreted, we use a maximum depth of
4.

## References

  * [1] Rosa Aghdam, Mojtaba Ganjali, Xiujun Zhang, and Changiz Eslahchi.  CN: a consensus algorithm for inferring gene regulatory networks using the sorder algorithm and conditional mutual information test.  Molecular BioSystems, 11(3):942–949, 2015.
  * [2] John Aitchison.  The statistical analysis of compositional data.  Journal of the Royal Statistical Society: Series B (Methodological), 44(2):139–160, 1982.
  * [3] Melis N. Anahtar, Jason H. Yang, and Sanjat Kanjilal.  Applications of machine learning to the problem of antimicrobial resistance: an emerging model for translational research.  Journal of Clinical Microbiology, 59(7):e01260–20, 2021.
  * [4] Michelle Badri, Zachary D Kurtz, Richard Bonneau, and Christian L Müller.  Shrinkage improves estimation of microbial associations under different normalization methods.  NAR genomics and bioinformatics, 2(4):lqaa100, 2020.
  * [5] Mariana Belgiu and Lucian Drăguţ.  Random forest in remote sensing: A review of applications and future directions.  ISPRS journal of photogrammetry and remote sensing, 114:24–31, 2016\.
  * [6] Roeland L. Berendsen, Corné M.J. Pieterse, and Peter A.H.M. Bakker.  The rhizosphere microbiome and plant health.  Trends in Plant Science, 17(8):478–486, 2012.
  * [7] Benjamin J Callahan, Paul J McMurdie, and Susan P Holmes.  Exact sequence variants should replace operational taxonomic units in marker-gene data analysis.  The ISME journal, 11(12):2639–2643, 2017.
  * [8] Tony F Chan, Gene H Golub, and Randall J LeVeque.  Algorithms for computing the sample variance: Analysis and recommendations.  The American Statistician, 37(3):242–247, 1983.
  * [9] Soumya Ghosh, Jiayu Yao, and Finale Doshi-Velez.  Model selection in bayesian neural networks via horseshoe priors.  Journal of Machine Learning Research, 20(182):1–46, 2019.
  * [10] Nicholas J Gotelli and Robert K Colwell.  Quantifying biodiversity: procedures and pitfalls in the measurement and comparison of species richness.  Ecology letters, 4(4):379–391, 2001.
  * [11] Ksenia Guseva, Sean Darcy, Eva Simon, Lauren V Alteio, Alicia Montesinos-Navarro, and Christina Kaiser.  From diversity to complexity: Microbial networks in soils.  Soil Biology and Biochemistry, 169:108604, 2022.
  * [12] José Miguel Hernández-Lobato and Ryan P. Adams.  Probabilistic backpropagation for scalable learning of Bayesian Neural Networks.  ICML’15: Proceedings of the 32nd International Conference on International Conference on Machine Learning, 37:1861–1869, July 2015.
  * [13] Athar Khodabakhsh, Tobias P. Loka, Sébastien Boutin, Dennis Nurjadi, and Bernhard Y. Renard.  Predicting decision-making time for diagnosis over ngs cycles: An interpretable machine learning approach.  bioRxiv, 2023.
  * [14] Zachary D Kurtz, Christian L Müller, Emily R Miraldi, Dan R Littman, Martin J Blaser, and Richard A Bonneau.  Sparse and compositionally robust inference of microbial ecological networks.  PLoS computational biology, 11(5):e1004226, 2015.
  * [15] Seongmin Lim, Jin-Hyung Kim, and Hae-Dong Kim.  Strategy for on-orbit space object classification using deep learning.  Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 235(15):2326–2341, 2021.
  * [16] Josep A Martín-Fernández, Carles Barceló-Vidal, and Vera Pawlowsky-Glahn.  Dealing with zeros and missing values in compositional data sets using nonparametric imputation.  Mathematical Geology, 35(3):253–278, 2003.
  * [17] Josep-Antoni Martín-Fernández, Karel Hron, Matthias Templ, Peter Filzmoser, and Javier Palarea-Albaladejo.  Bayesian-multiplicative treatment of count zeros in compositional data sets.  Statistical Modelling, 15(2):134–158, 2015.
  * [18] Paul J McMurdie and Susan Holmes.  Waste not, want not: why rarefying microbiome data is inadmissible.  PLoS computational biology, 10(4):e1003531, 2014.
  * [19] Radford M. Neal.  MCMC Using Hamiltonian Dynamics.  CRC Press, 2011.
  * [20] Radford M. Neal.  Bayesian learning for neural networks, volume 118.  Springer Science & Business Media, 2012.
  * [21] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [22] Mariana Buongermino Pereira, Mikael Wallroth, Viktor Jonsson, and Erik Kristiansson.  Comparison of normalization methods for the analysis of metagenomic gene abundance data.  BMC Genomics, 19(1):274, 2018.
  * [23] Stefanie Peschel, Christian L Müller, Erika von Mutius, Anne-Laure Boulesteix, and Martin Depner.  Netcomi: network construction and comparison for microbiome data in r.  Briefings in bioinformatics, 22(4):bbaa290, 2021.
  * [24] JF Mandujano Reyes, E Walleser, S Hachenberg, S Gruber, M Kammer, C Baumgartner, R Mansfeld, K Anklam, and D Döpfer.  Full model selection using regression trees for numeric predictions of biomarkers for metabolic challenges in dairy cows.  Preventive Veterinary Medicine, 193:105422, 2021.
  * [25] Quan Sun, Bernhard Pfahringer, and Michael Mayo.  Full model selection in the space of data mining operators.  In Proceedings of the 14th annual conference companion on genetic and evolutionary computation, pages 1503–1504, 2012.
  * [26] Quan Sun, Bernhard Pfahringer, and Michael Mayo.  Towards a framework for designing full model selection and optimization systems.  In International Workshop on Multiple Classifier Systems, pages 259–270. Springer, 2013.
  * [27] Pankaj Trivedi, Manuel Delgado-Baquerizo, Ian C. Anderson, and Brajesh K. Singh.  Response of soil properties and microbial communities to agriculture: Implications for primary productivity and soil health indicators.  Frontiers in Plant Science, 7, 2016.
  * [28] Cameron Wagg, Klaus Schlaeppi, Samiran Banerjee, Eiko E Kuramae, and Marcel GA van der Heijden.  Fungal-bacterial diversity and microbiome complexity predict ecosystem functioning.  Nature communications, 10(1):4841, 2019.
  * [29] Yinglin Xia, Jun Sun, Ding-Geng Chen, et al.  Statistical analysis of microbiome data with R, volume 847.  Springer, 2018.
  * [30] Tao Yan, Shui-Long Shen, Annan Zhou, and Xiangsheng Chen.  Prediction of geological characteristics from shield operational parameters by integrating grid search and k-fold cross validation into stacking classification algorithm.  Journal of Rock Mechanics and Geotechnical Engineering, 2022.
  * [31] Jiayu Yao, Weiwei Pan, Soumya Ghosh, and Finale Doshi-Velez.  Quality of uncertainty quantification for bayesian neural network inference.  Proceedings at the International Conference on Machine Learning: Workshop on Uncertainty & Robustness in Deep Learning, June 2019.
  * [32] In-Kwon Yeo and Richard A Johnson.  A new family of power transformations to improve normality or symmetry.  Biometrika, 87(4):954–959, 2000.
  * [33] Grace Yoon, Irina Gaynanova, and Christian L Müller.  Microbial networks in spring-semi-parametric rank-based correlation and partial correlation estimation for quantitative microbiome data.  Frontiers in genetics, 10:516, 2019.

## Data availability

Data are available upon request to Richard Lankau (lankau@wisc.edu).

## Code Availability

All reproducible scripts are open source and publicly available in
https://github.com/solislemuslab/soil-microbiome-nn.

## Acknowledgements

The authors thank Linda Kinkel and the Kinkel group. This work was supported
by the Department of Energy [DE-SC0021016 to CSL]. The work was also supported
by USDA Specialty Crop Multi-State Grant Program award SCMP1701.

## Author contributions

CSL and RL developed the idea. RL and SS collected the data. XT and RL led all
statistical analyses from data preprocessing to fitting of machine learning
models, as well as summarizing the results by the creation of figures. CSL, XT
and RA wrote the initial complete draft of the manuscript. RL and SS
contributed in biological interpretations and edition of the manuscript. All
authors read and approved the final manuscript.

## Competing interests

The authors declare that they have no competing interests.

Supplementary Material:

Human limits in Machine Learning: Prediction of plant phenotypes using soil
microbiome data

Soil health |
---|---
Physicochemical (Soil) | pH
| Texture
| Organic matter content
| Cation/Anion exchange capacity
| Carbon fractions
| Organic nitrogen
| Available macronutrients
| Available micronutrients
| Potential carbon mineralization rate
Microbiome (OTU and alpha diversity) | Bacterial community composition, structure, diversity
| Fungal community composition, structure, diversity
| Bacterial and fungal population abundance
Pathogen | Verticillium dahliae, Pathogenic Streptomyces population
| Disease suppressiveness (DS)
| lab-estimated soil disease suppressiveness
Potato crop assessment |
Yield | Tuber Yield, yield by size class
Disease severity |
| Verticillium dahlia incidence
| Common scab disease incidence and severity
| Other diseases including hollow heart, silver scurf and black
| scurf
Management practices |
Field history | Crop rotation
Practices around growing season | Potato cultivar selection, fumigation,
| cover crop, fertilization, pesticide, irrigation
Table 1: Description of variables included in the soil dataset. Input: All
samples per varietySubset of samples only for Russet varietyCompute the median
of the yield valuesFor each sample in Russet varietyYield value $>$
medianLabel = 1Label = 0yesno Figure 1: Flowchart for binarizing the
continuous yield response into binary labels for the Russet variety. The same
procedure is used for every variety in the dataset. Input: training set of
original OTUs & all labels# Label $0$ $<$ 400Select Label $0$# Label $1$ $<$
400Select Label $1$For each potato varietySelect the subset of samples for
that varietyCalculate the mean and variance of each OTU in the subsetGenerate
noises $\mathcal{N}\sim(\mu/100,\sigma/100)$Augmented value = max(0,
original+noise)Add this augmented sampleOutput: 800 samplesyesnoyesno Figure
2: Flowchart with the data augmentation algorithm. The target sample size for
the training set is 800 with 400 samples for each label. The noise that we
artificially generate needs to be variety-specific before adding to the
original samples so that the biological implications of the original samples
would be preserved. Figure 3: Boxplots of weighted F1 scores for random forest
models for different types of predictors (rows) and different yield or disease
outcomes (columns). For description of the rows, see Table 3 in the main text.
The range of each boxplot depicts the weighted F1 scores for datasets in
different taxonomic levels and different normalization and zero replacement
strategies. The dashed line corresponds to the weighted F1 score when fitting
the model with random datasets(see Section B.2.1). Figure 4: Boxplots of
weighted F1 scores for Bayesian neural network models for different types of
predictors (rows) and different yield or disease outcomes (columns). For
description of the rows, see Table 3 in the main text. The range of each
boxplot depicts the weighted F1 scores for datasets in different taxonomic
levels and different normalization and zero replacement strategies. The dashed
line corresponds to the weighted F1 score when fitting the model with random
datasets (see Section B.2.1).

| KBest | Mutual | LR | DT | GB | RF | Max | TOTAL
---|---|---|---|---|---|---|---|---
Firmicutes | X | X | X | X | X | X | X | 7
Patescibacteria | X | X | X | X | X | X | X | 7
Myxococcota | - | - | X | X | X | X | X | 5
Methylomirabilota | X | X | X | X | - | - | X | 5
Acidobacteriota | - | - | X | X | X | X | X | 5
Verrucomicrobiota | - | - | X | X | X | - | X | 4
Chloroflexi | - | - | X | X | - | X | X | 4
Desulfobacterota | X | X | X | - | - | X | - | 4
NB1.j | X | X | - | X | - | - | - | 3
Thermoplasmatota | X | X | - | - | X | - | - | 3
RCP2.54 | - | - | - | X | X | X | - | 3
WS2 | - | X | - | X | X | - | - | 3
Dependentiae | - | - | - | X | X | X | - | 3
WPS.2 | X | X | - | - | - | X | - | 3
Cyanobacteria | - | - | X | - | X | - | X | 3
Armatimonadota | X | X | - | X | - | - | - | 3
Elusimicrobiota | X | - | - | - | - | X | - | 2
Latescibacterota | X | - | - | - | - | X | - | 2
Bdellovibrionota | - | - | - | - | X | X | - | 2
Bacteroidota | - | - | X | - | - | - | X | 2
Spirochaetota | X | X | - | - | - | - | - | 2
Nitrospirota | - | - | X | - | - | - | X | 2
Actinobacteriota | - | - | X | - | - | - | X | 2
Planctomycetota | - | - | - | X | - | X | - | 2
Proteobacteria | - | - | X | - | - | - | X | 2
SAR324 | - | X | - | - | - | - | - | 1
Abditibacteriota | - | - | - | - | X | - | - | 1
Nitrospinota | - | X | - | - | - | - | - | 1
Nanoarchaeota | - | X | - | - | - | - | - | 1
Gemmatimonadota | - | - | - | - | - | - | X | 1
GAL15 | - | - | - | - | X | - | - | 1
Deinococcota | X | - | - | - | - | - | - | 1
MBNT15 | X | - | - | - | - | - | - | 1

Table 2: The Operational Taxonomic Units (OTUs) based on their maximum values
in samples are sorted and the features that their maximum values are among the
top 30% (Max column) are selected. The selected features by the SelectKBest
method are marked in the KBest column. Columns three to six are the result of
applying of logistic regression (LR), decision tree (DT), Gradient Boosting
(GB), or Random Forrest (RF) as the choice of the algorithm using in the
Recursive Feature Elimination method. The selected features are also
considered based on mutual information statistics as shown in the Mutual
column. We assigned a TOTAL value to each OTU based on the number of times the
OTU gets picked by any of the seven criteria (six ML models and the max OTU
value) shown in the TOTAL column for the Phylum level associated with Scabpit
disease. The OTUs with scoring values higher than zero are shown.

| Degree 0 | Degree 1 | Degree difference
---|---|---|---
MBNT15 | 1 | 9 | 8
Actinobacteriota | 8 | 3 | 5
NB1.j | 3 | 8 | 5
Sva0485 | 6 | 2 | 4
Bacteroidota | 4 | 7 | 3
Halobacterota | 4 | 1 | 3
Methylomirabilota | 5 | 2 | 3
Proteobacteria | 3 | 6 | 3
Chloroflexi | 4 | 6 | 2
Cyanobacteria | 0 | 2 | 2
Entotheonellaeota | 3 | 1 | 2
GAL15 | 1 | 3 | 2
Hydrogenedentes | 2 | 0 | 2
Patescibacteria | 0 | 2 | 2
Planctomycetota | 0 | 2 | 2
RCP2.54 | 2 | 4 | 2
WPS.2 | 2 | 0 | 2
Zixibacteria | 1 | 3 | 2
Abditibacteriota | 0 | 1 | 1
Desulfobacterota | 4 | 5 | 1
Elusimicrobiota | 0 | 1 | 1
Fibrobacterota | 2 | 3 | 1
Gemmatimonadota | 2 | 1 | 1
Myxococcota | 3 | 2 | 1
Nanoarchaeota | 2 | 1 | 1
Nitrospirota | 3 | 2 | 1
Spirochaetota | 3 | 4 | 1
Sumerlaeota | 1 | 0 | 1
Thermoplasmatota | 1 | 0 | 1
Verrucomicrobiota | 0 | 1 | 1
Acidobacteriota | 3 | 3 | 0
Armatimonadota | 0 | 0 | 0
Bdellovibrionota | 2 | 2 | 0
Crenarchaeota | 1 | 1 | 0
Deinococcota | 0 | 0 | 0
Dependentiae | 0 | 0 | 0
FCPU426 | 0 | 0 | 0
Firmicutes | 1 | 1 | 0
Latescibacterota | 2 | 2 | 0
Nitrospinota | 4 | 4 | 0
SAR324 | 0 | 0 | 0
WS2 | 1 | 1 | 0

Table 3: List of OTUs in microbial networks constructed by SPRING for two classes of scabpit. Degree $0$ correspond to the degree of the OTU node in the network built with the samples with Label $0$. Degree $1$ corresponds to the degree of the OTU node in the network built with samples with Label $1$ Parameters | Description | Tuning Parameters
---|---|---
n_estimators | The number of trees in the forest. | [100,200,500]
min_samples_split | Minimum necessary number of samples to split an internal node | [8,10]
min_samples_leaf | The minimum needed the number of samples at a leaf node. | [3,4,5]
max_depth | The maximum depth of the tree | [80,90]
criterion | The function used to evaluate a split’s quality | (’gini’,’entropy’)
Table 4: Different values for parameters for Random Forest Model. Figure 5:
Scatterplot for F1 scores in the non-diseased (Label $0$) group versus the F1
scores in the diseased (Label $1$) group by disease outcome (rows) and
randomization strategy (columns). Each panel has 200 generated datasets (each
blue point is a random dataset), and in a red point, we should the results on
the real data. These results are for data at the Phylum level. Fewer points in
scab and black scurf diseases are due to extreme imbalance. For example, there
are only 11 samples without disease (Label $0$) out of 46 testing samples in
Scab which results in 11 label 0 samples and only 12 possible F1 values for
both 0 and 1 classes, i.e, only 12 points. Since there are far more label 1
samples in Scab, the model is biased towards label 1 and tend to predict most
testing samples as label 1, so the F1 values are clustered around 75-85% for
label 1. In this case, the points are more likely to stack on top of one
another, thus leading to fewer points in the visualization.  Figure 6:
Weighted F1 scores (y-axis) for Random Forest and Bayesian Neural Network (NN)
models for yield by plant under the 20 normalization/zero replacement
strategies (x-axis). The lack of pattern prevents us from making
recommendations of optimal strategies for microbiome OTU data. We can
conclude, however, that taxonomic levels, normalization and zero replacement
strategies have an effect on the prediction accuracy of the models as
evidenced by the broad range displayed by the points. Figure 7: Weighted F1
scores (y-axis) for Random Forest and Bayesian Neural Network (NN) models for
yield by meter under the 20 normalization/zero replacement strategies
(x-axis). The lack of pattern prevents us from making recommendations of
optimal strategies for microbiome OTU data. We can conclude, however, that
taxonomic levels, normalization and zero replacement strategies have an effect
on the prediction accuracy of the models as evidenced by the broad range
displayed by the points. Figure 8: Weighted F1 scores (y-axis) for Random
Forest and Bayesian Neural Network (NN) models for black scurf disease under
the 20 normalization/zero replacement strategies (x-axis). The lack of pattern
prevents us from making recommendations of optimal strategies for microbiome
OTU data. We can conclude, however, that taxonomic levels, normalization and
zero replacement strategies have an effect on the prediction accuracy of the
models as evidenced by the broad range displayed by the points. Figure 9:
Weighted F1 scores (y-axis) for Random Forest and Bayesian Neural Network (NN)
models for scab disease under the 20 normalization/zero replacement strategies
(x-axis). The lack of pattern prevents us from making recommendations of
optimal strategies for microbiome OTU data. We can conclude, however, that
taxonomic levels, normalization and zero replacement strategies have an effect
on the prediction accuracy of the models as evidenced by the broad range
displayed by the points. Figure 10: Weighted F1 scores (y-axis) for Random
Forest and Bayesian Neural Network (NN) models for superficial scab
(Scabsuper) disease under the 20 normalization/zero replacement strategies
(x-axis). The lack of pattern prevents us from making recommendations of
optimal strategies for microbiome OTU data. We can conclude, however, that
taxonomic levels, normalization and zero replacement strategies have an effect
on the prediction accuracy of the models as evidenced by the broad range
displayed by the points. Figure 11: Weighted F1 scores (y-axis) by Random
Forest and Bayesian Neural Network (NN) models for yield by plant by feature
selection strategy (x-axis) including all OTUs (All-OTU), OTUs selected by the
ML method (OTU-S1), the network comparison method (OTU-S2), both methods
(OTU-S3), or neither method (OTU-S0). Figure 12: Weighted F1 scores (y-axis)
by Random Forest and Bayesian Neural Network (NN) models for yield by meter by
feature selection strategy (x-axis) including all OTUs (All-OTU), OTUs
selected by the ML method (OTU-S1), the network comparison method (OTU-S2),
both methods (OTU-S3), or neither method (OTU-S0). Figure 13: Weighted F1
scores (y-axis) by Random Forest and Bayesian Neural Network (NN) models for
black scurf disease by feature selection strategy (x-axis) including all OTUs
(All-OTU), OTUs selected by the ML method (OTU-S1), the network comparison
method (OTU-S2), both methods (OTU-S3), or neither method (OTU-S0). Figure 14:
Weighted F1 scores (y-axis) by Random Forest and Bayesian Neural Network (NN)
models for scab disease by feature selection strategy (x-axis) including all
OTUs (All-OTU), OTUs selected by the ML method (OTU-S1), the network
comparison method (OTU-S2), both methods (OTU-S3), or neither method (OTU-S0).
Figure 15: Weighted F1 scores (y-axis) by Random Forest and Bayesian Neural
Network (NN) models for superficial scab (Scabsuper) disease by feature
selection strategy (x-axis) including all OTUs (All-OTU), OTUs selected by the
ML method (OTU-S1), the network comparison method (OTU-S2), both methods
(OTU-S3), or neither method (OTU-S0). Figure 16: Full model selection decision
tree with a maximum depth of 4 summarizing the results of Random Forest models
on yield by plant. When the condition at a node is true, we follow the branch
on the left, and when the condition is false, we follow the branch on the
right. The percentage of data preprocessing options and mean of weighted F1
scores are shown in each nodes. Figure 17: Full model selection decision tree
with a maximum depth of 4 summarizing the results of Random Forest models on
yield by meter. When the condition at a node is true, we follow the branch on
the left, and when the condition is false, we follow the branch on the right.
The percentage of data preprocessing options and mean of weighted F1 scores
are shown in each nodes. Figure 18: Full model selection decision tree with a
maximum depth of 4 summarizing the results of Random Forest models on black
scurf disease. When the condition at a node is true, we follow the branch on
the left, and when the condition is false, we follow the branch on the right.
The percentage of data preprocessing options and mean of weighted F1 scores
are shown in each nodes. Figure 19: Full model selection decision tree with a
maximum depth of 4 summarizing the results of Random Forest models on scab
disease. When the condition at a node is true, we follow the branch on the
left, and when the condition is false, we follow the branch on the right. The
percentage of data preprocessing options and mean of weighted F1 scores are
shown in each nodes. Figure 20: Full model selection decision tree with a
maximum depth of 4 summarizing the results of Random Forest models on
superficial scab (Scabsuper) disease. When the condition at a node is true, we
follow the branch on the left, and when the condition is false, we follow the
branch on the right. The percentage of data preprocessing options and mean of
weighted F1 scores are shown in each nodes. Figure 21: Full model selection
decision tree with a maximum depth of 4 summarized the results of Bayesian
Neural Network models on yield by plant. When the condition at a node is true,
we follow the branch on the left, and when the condition is false, we follow
the branch on the right. The percentage of data preprocessing options and mean
of weighted F1 scores are shown in each nodes. Figure 22: Full model selection
decision tree with a maximum depth of 4 summarized the results of Bayesian
Neural Network models on yield by meter. When the condition at a node is true,
we follow the branch on the left, and when the condition is false, we follow
the branch on the right. The percentage of data preprocessing options and mean
of weighted F1 scores are shown in each nodes. Figure 23: Full model selection
decision tree with a maximum depth of 4 summarized the results of Bayesian
Neural Network models on black scurf disease. When the condition at a node is
true, we follow the branch on the left, and when the condition is false, we
follow the branch on the right. The percentage of data preprocessing options
and mean of weighted F1 scores are shown in each nodes. Figure 24: Full model
selection decision tree with a maximum depth of 4 summarized the results of
Bayesian Neural Network models on scab disease. When the condition at a node
is true, we follow the branch on the left, and when the condition is false, we
follow the branch on the right. The percentage of data preprocessing options
and mean of weighted F1 scores are shown in each nodes. Figure 25: Full model
selection decision tree with a maximum depth of 4 summarized the results of
Bayesian Neural Network models on superficial scab (Scabsuper) disease. When
the condition at a node is true, we follow the branch on the left, and when
the condition is false, we follow the branch on the right. The percentage of
data preprocessing options and mean of weighted F1 scores are shown in each
nodes. Figure 26: Boxplots of the weighted F1 scores by Random Forest (RF) and
Bayesian Neural Network (BNN) models for yield by plant using environmental
predictors. Figure 27: Boxplots of the weighted F1 scores by Random Forest
(RF) and Bayesian Neural Network (BNN) models for yield by meter using
environmental predictors. Figure 28: Boxplots of the weighted F1 scores by
Random Forest (RF) and Bayesian Neural Network (BNN) models for black scurf
disease using environmental predictors. Figure 29: Boxplots of the weighted F1
scores by Random Forest (RF) and Bayesian Neural Network (BNN) models for scab
disease outcome using environmental predictors. Figure 30: Boxplots of the
weighted F1 scores by Random Forest (RF) and Bayesian Neural Network (BNN)
models for superficial scab (Scabsuper) disease using environmental
predictors. Figure 31: Boxplots with the weighted F1 scores (y-axis) by Random
Forest and Bayesian Neural Network (NN) models (x-axis) for yield per plant.
The models including both types of predictors outperform other models, yet
models including OTU data alone (OTU-S3) are comparable which suggests that
the microbial information indeed contains signal to predict the disease
outcome on its own. However, OTU data is more expensive to collect, and
perhaps not necessary, given that the model without OTU data (Soil) performs
just as accurately (blue dashed line). Figure 32: Boxplots with the weighted
F1 scores (y-axis) by Random Forest and Bayesian Neural Network (NN) models
(x-axis) for yield per meter. The models including both types of predictors
outperform other models, yet models including OTU data alone (OTU-S3) are
comparable which suggests that the microbial information indeed contains
signal to predict the disease outcome on its own. However, OTU data is more
expensive to collect, and perhaps not necessary, given that the model without
OTU data (Soil) performs just as accurately (blue dashed line). Figure 33:
Boxplots with the weighted F1 scores (y-axis) by Random Forest and Bayesian
Neural Network (NN) models (x-axis) for black scurf disease. The models
including both types of predictors outperform other models, yet models
including OTU data alone (OTU-S3) are comparable which suggests that the
microbial information indeed contains signal to predict the disease outcome on
its own. However, OTU data is more expensive to collect, and perhaps not
necessary, given that the model without OTU data (Soil) performs just as
accurately (blue dashed line). Figure 34: Boxplots with the weighted F1 scores
(y-axis) by Random Forest and Bayesian Neural Network (NN) models (x-axis) for
scab disease. The models including both types of predictors outperform other
models, yet models including OTU data alone (OTU-S3) are comparable which
suggests that the microbial information indeed contains signal to predict the
disease outcome on its own. However, OTU data is more expensive to collect,
and perhaps not necessary, given that the model without OTU data (Soil)
performs just as accurately (blue dashed line). Figure 35: Boxplots with the
weighted F1 scores (y-axis) by Random Forest and Bayesian Neural Network (NN)
models (x-axis) for superficial scab (Scabsuper) disease. The models including
both types of predictors outperform other models, yet models including OTU
data alone (OTU-S3) are comparable which suggests that the microbial
information indeed contains signal to predict the disease outcome on its own.
However, OTU data is more expensive to collect, and perhaps not necessary,
given that the model without OTU data (Soil) performs just as accurately (blue
dashed line).