Figure 1. Histograms and density distributions of the averaged (a) atomic masses and (b) electronegativities based on the chemical compositions in the widely used datasets in MI.

2.1 Data preparation and construction of the machine learning model

We used experimental property data from Starrydata2, a material database developed by our group [,], comprising 42,005 samples and 2,236,338 records extracted from 7,698 papers as on 27 September 2021 []. This is the largest dataset that includes the temperature dependence of the experimental properties in the field of thermoelectric materials, worldwide. Although there is a data bias toward the thermoelectric field (), this dataset has a comprehensive collection of experimental data for various material families, performances, and temperature ranges, which reduces the selection bias compared to other experimental value databases [,]. All the data from Starrydata2 are available for free on GitHub []. This includes two types of data: The first is raw data extracted directly from the figures in a paper; the second includes data extracted only from the physical property values where the x-axis denotes the temperature, and each property value is interpolated by a 5-th order polynomial for every 50 K temperature point. In this study, the data in Starrydata2 as on 27 September 2021, interpolated by the 5-th order polynomial, were used for training and validation. In particular, only thermoelectric property data (294,616 records) from 6,594 papers on thermoelectric materials were used in this study.

The performance of thermoelectric materials was evaluated as a dimensionless figure of merit, zT=S²σ/κ, where S is the Seebeck coefficient, σ is the electrical conductivity, κ is the thermal conductivity, and T is the temperature. The discovery of high-performance thermoelectric materials involves searching for materials with a high-power factor (PF = S²σ) and low κ. For several years, degenerate semiconductors based on heavy metals such as Bi, Pb, Te, and Sb have been the mainstream materials exhibiting high zT. Numerous other material systems, such as skutterudite and clathrate compounds, Zintl phase and Heusler compounds, and oxides and silicide compounds, have been reported as thermoelectric materials and are included in the dataset used in this study.

Because the data in Starrydata2 are extracted from previously published papers, numerous errors exist in the papers themselves and in the extraction process. Therefore, the data were strictly preprocessed using the following two procedures. These procedures are expected to provide reasonable values for T, S, σ, κ, and zT.

1. Remove the records in which at least T, S, σ, κ, or zT is missing (63,405 records).
2. Remove the records for which the mean absolute percentage error rate (MAPE) between zT_calc calculated from T, S, σ, and κ, and zT extracted directly from the paper is greater than 5% (35,393 records).

The input vectors used for ML included 26 feature values based on the chemical composition and the (measured) temperature T. The features were the calculated mean, variance, and difference (maximum value-minimum value within the major elements) of the ‘group’, ‘period’, ‘atomic number’, ‘Mendeleev number’, ‘atomic radius’, ‘atomic weight’, ‘electro negativity’, and ‘VEC’ of the containing elements, and ‘number of containing elements’ and ‘number of major elements’. Here, the major elements were those with a ratio of 0.1 or greater in the chemical composition such that the overall ratio was unity []. It is interesting to note that as Starrydata2 records the temperature dependence of the experimental physical properties, the temperature can be added to the input vector. Features other than the temperature, are similar to the Magpie/Matminer [,] feature set. However, for the material in Starrydata2, several small amounts of additive elements are often introduced to modulate the carrier concentration or reduce lattice thermal conductivity. Moreover, there are complex chemical compositions consisting of up to 10 elements. The ‘range’, ‘minimum’, and ‘maximum’ attributes in the Magpie/Matminer feature set may overestimate the effect of these additive elements, and were therefore excluded from the input vector. In addition, the electronic structure attributes and ab initio calculated attributes were excluded from our feature set to avoid overcomplicating the feature space. Records that could not be converted into feature vectors were deleted (35,332 records). The Python library for material analysis, pymatgen [], was used to create records of various features. The target variables were thermoelectric properties S, σ, κ, and zT_calc. The training and test data were split based on the publication year of the paper rather than the usual random split. Data published before 2020 were used as training data (21,775 records), whereas those from 2021 onward were used as test data (563 records). Thus, the problem was set as predicting materials in 2021 that were unknown before 2020. Because the materials up to 2020, which would normally have been included in 2021, are completely removed, the test data include only completely unknown chemical compositions. Note that the test data is used to evaluate the generalization performance, and is different from the validation data used to optimize the hyperparameters.

For the ML method, we used the random forest algorithm. This is an ensemble learning algorithm, which combines multiple decision trees to improve the generalization performance. In addition, it is robust to noise and overfitting, and has often been used in previous research in materials science [,]. The hyperparameters were optimized using the grid search provided by Scikit-learn []. The MAPE and root mean squared logarithmic error (RMSLE) were used to evaluate the prediction accuracy.(1)MAPE=100n∑i=1nti−yiyi,(1)(2)RMSLE=1n∑i=1nlogti+1−logyi+12,(2)

where n is the number of test data points, t is the experimental value, and y is the predicted value. As target variables S and σ are physical properties with a wide range of variation in the order of magnitude, the absolute error and root mean square error, which are generally used, are not used in this study.

In order to propose new materials, the ML model constructed in this study was applied to the chemical compositions recorded in the Materials Project [] to predict various thermoelectric properties. We used all the Materials Project data as on 18 October 2018 (83,989 records). Note that only the chemical compositions of stable materials (energy above hull > = 15 meV), including known metastable materials, and the materials included in the AD were used.

2.2 Clustering of material families

To clarify the material space of the dataset, soft clustering using variational Bayesian estimation of Gaussian mixture distribution was used. In materials science, the use of soft clustering is reasonable because some materials belong to multiple material family clusters. For variational Bayesian estimation of the mixture Gaussian distribution, we used an algorithm provided by the Python ML library, Scikit-learn []. The variance-covariance matrix of the Gaussian distribution was used as the common variance-covariance matrix for each cluster (tied type) [], and the number of clusters was set to 15. The input vector used for clustering is the vector described in . The number of clusters was arbitrarily determined as a value that could be identified by thermoelectric material experts. To clarify the characteristics of each cluster, the elements contained in the materials within the clusters were counted, and the top three elements with the highest number of occurrences were assigned labels.

2.3 Applicability domain using the k-nearest neighbor method

There are various types of ADs, such as Bounding Box, Convex Hull, Clustering, and the k-nearest neighbor method, depending on how the domain is defined [,]. Bounding Box is the simplest and fastest method for defining an n-dimensional hyper-rectangle as an AD based on the maximum and minimum values of each descriptor. However, due to the simple box of the defined feature space, several false positive domains are generated when the data points are nonuniformly distributed. Convex Hull is a method that improves on the Bounding Box using the smallest convex domain of the feature space as the AD; however, false positives may occur in concave regions when non-convex domains exist in the feature space. In addition, this method is computationally expensive. AD using Clustering, such as k-means, is computationally less expensive, and the definition of the AD is not as crude as that of the Bounding Box. Therefore, this method is effective for defining the AD in cases where reasonable clustering can be achieved. Moreover, in terms of clearly defining the learned range by Clustering, it is similar to the leave-one-cluster-out (LOCO) cross-validation proposed by Bryce Meredig et al []. Therefore, it is also an effective tool for evaluating extrapolation. However, false positive domains may occur when the training data points are non-uniformly distributed and the cluster center of the gravity point is not in a data dense domain. AD with k-nearest neighbors is a method that defines regions using the average distance from the k nearest neighbors in all training data points as the threshold. This method is more rigorous than its other counterparts because the AD is clearly defined by the distance between each training data and the test data. In this study, we focused on the AD with the k-nearest neighbor method to rigorously discuss the effects of data bias. Clustering in was not used to define AD but only to identify material families for the purpose of providing human interpretability.

Among the k-nearest neighbor methods, we used the approach proposed by Sahigara et al []. In general, a single threshold (e.g. the 95th percentile) is determined based on the average distance around k for all the data points. However, a single threshold may be determined to be biased toward data points in regions of high data point density (e.g. existing high-performance materials) and may not fit regions of low data point density (e.g. new materials). In the methods proposed by Sahigara et al., thresholds are set for each data point, and AD is determined individually. In addition, the value of the number of neighbors, k, is important for determining the AD using the k-neighborhood method. A low value of k leads to overly strict limitations, while a high value unnecessarily expands the AD. Therefore, k is determined by maximizing the percentage of validation data retained in the AD using the Monte Carlo method. In this study, we determined the optimal k by randomly dividing the AD 1,000 times such that 20% (4,355 records) of the training data (21,775 records) became the validation data, and obtained the percentage of data retained in the AD for each different k.

Figure 2. Dependence of the number of training data on the publication year of the paper. The number of training data is stacked with the results of material family clustering.

Figure 3. Box-And-Whisker plot of the test samples (%) retained within the AD for different k-values during k-optimization.

Figure 4. Actual vs. predicted TE properties for the test data. The darker the plot color, the more are the number of neighboring known materials in the AD.

Figure 5. Dependence of (a) MAPE (line chart) and (b) RMSLE (line chart) on the number of adjacent known materials in the AD and the number of test data (unknown materials) on the stacked surface graph (second axis). The colors of the faces are identical to those of the material system clusters in .

Figure 6. Top-10 predicted S values at each temperature (300 K–900 K) for the chemical composition of the Materials Project data included in the AD. The chemical compositions contained in the training data are masked in gray. At the bottom of each block, the cluster to which the chemical composition belongs, the number of adjacent known materials, and the range of the predicted S calculated from the MAPE are shown.