Identification of taxon through classification with partial reject options

TL;DR

This paper introduces a Bayesian classification method for identifying taxa using mixed trait data, allowing for uncertainty and partial rejection, with applications demonstrated on ornithological datasets.

Contribution

It presents a novel Bayesian discriminant analysis with partial reject options and outlier safeguarding, applicable to mixed continuous and ordinal traits.

Findings

Effective classification with uncertainty handling

Outlier detection using Bayesian p-value analogue

Validated on original ornithological datasets

Abstract

Identification of taxa can significantly be assisted by statistical classification based on trait measurements in two major ways; either individually or by phylogenetic (clustering) methods. In this paper we present a general Bayesian approach for classifying species individually based on measurements of a mixture of continuous and ordinal traits as well as any type of covariates. It is assumed that the trait vector is derived from a latent variable with a multivariate Gaussian distribution. Decision rules based on supervised learning are presented that estimate model parameters through blockwise Gibbs sampling. These decision regions allow for uncertainty (partial rejection), so that not necessarily one specific category (taxon) is output when new subjects are classified, but rather a set of categories including the most probable taxa. This type of discriminant analysis employs reward functions with a set-valued input argument, so that an optimal Bayes classifier can be defined. We also present a way of safeguarding against outlying new observations, using an analogue of a $p$-value within our Bayesian setting. Our method is illustrated on an original ornithological data set of birds. We also incorporate model selection through cross-validation, examplified on another original data set of birds.

Tables

Table 1. Table 1: The distribution of the Acrocephalus data over species, covariates and traits. Note that the Paddyfield Warbler has few observations overall, but almost all of them are complete, whereas the Reed Warbler has many observations but only a small percentage of them are complete.

		Species
Trait	Covariate	Reed W.	Blyth’s Reed W.	Paddyfield W.	Marsh W.
Wing	Juvenile	36 023	41	19	2 534
length	Adult	14 536	68	12	870
Notch	Juvenile	830	41	18	2 212
length	Adult	472	68	12	579
Notch	Juvenile	410	41	19	416
Position	Adult	77	68	12	40

Table 2. Table 2: Hyperparameter values for the prior of the regression parameter matrices. These values are informed by Malmhagen et al. ( 2013 ) , except for the first element of the second row of 𝐁 30 subscript 𝐁 30 \mathbf{B}_{30} and 𝐁 40 subscript 𝐁 40 \mathbf{B}_{40} , which we put to 0.75 0.75 0.75 , since we strongly believe that the pattern of slightly longer wings in adult plumages also hold for these two species.

		Reed Warbler		Marsh Warbler
	${𝐁}_{10}$		${𝐁}_{20}$
		Paddyfield Warbler		Blyth’s Reed Warbler
	${𝐁}_{30}$		${𝐁}_{40}$

Table 3. Table 3: Bayes estimates 𝐁 ^ i ​ m ​ k = 𝔼 ^ ​ ( 𝐁 i ​ m ​ k | 𝒟 i ) subscript ^ 𝐁 𝑖 𝑚 𝑘 ^ 𝔼 conditional subscript 𝐁 𝑖 𝑚 𝑘 subscript 𝒟 𝑖 \hat{\mathbf{B}}_{imk}=\hat{\mathbb{E}}(\mathbf{B}_{imk}|{\cal D}_{i}) and 𝚺 ^ i ​ k ​ l α = 𝔼 ^ ​ ( 𝚺 i ​ k ​ l α | 𝒟 i ) superscript subscript ^ 𝚺 𝑖 𝑘 𝑙 𝛼 ^ 𝔼 conditional superscript subscript 𝚺 𝑖 𝑘 𝑙 𝛼 subscript 𝒟 𝑖 \hat{\mathbf{\Sigma}}_{ikl}^{\alpha}=\hat{\mathbb{E}}(\mathbf{\Sigma}_{ikl}^{\alpha}|{\cal D}_{i}) of all parameters of the Acrocephalus model (cf. ( 7 )), rounded to two decimals, except the notch position trait where we present the categories that the regression parameters estimates fall into. The coding of these categories is explained in Table 2 of Malmhagen et al. ( 2013 ) .

		Reed warbler		Marsh warbler
	${\widehat{𝐁}}_1$		${\widehat{𝐁}}_2$
		Paddyfield Warbler		Blyth’s Reed Warbler
	${\widehat{𝐁}}_3$		${\widehat{𝐁}}_4$

Table 4. Table 4: The misclassification (as in rejecting the correct category) rate for each category using the uniform prior and reward function R ~ ~ 𝑅 \tilde{R} . The majority of observations are partially observed leading to the high rates of misclassification (cf. Section 4 ).

Species	$\rho =1$	$\rho =0.3$
Reed warbler	17.4%	5.1%
Blyth’s reed warbler	2.8%	0.9%
Paddyfield warbler	12.9%	3.2%
Marsh warbler	6.9%	2.2%

Table 5. Table 5: Misclassification rate is the probability of classifying a new observation wrongly, estimated by 100-fold cross validation, when new observations are assumed to have a uniform distribution over species.

Misclassification rate	Method
Age class	mclust	K-H
Adult	3.05 %	3.56 %
Juvenile	1.59 %	3.74 %

Table 6. Table 6: Classification table for adult birds using the cvMclustDA cross-validation method (referred to as mclust below), and the method proposed in this paper, (referred to as K-H below). The similarities of the methods are high, with a slightly higher accuracy for the mclust -method for Blyth’s reed warbler, and a slightly higher accuracy for the proposed method for the Paddyfield warbler.

	True species
	Reed warbler		Blyth’s reed w.		Paddyfield w.		Marsh warbler
Predicted species	mclust	K-H	mclust	K-H	mclust	K-H	mclust	K-H
Reed warbler	76	76	0	0	0	0	0	0
Blyth’s reed warbler	0	1	65	63	2	1	0	0
Paddyfield warbler	0	0	3	5	10	11	0	0
Marsh warbler	1	0	0	0	0	0	40	40

Table 7. Table 7: Classification table for juveline birds using the cvMclustDA cross-validation method (referred to as mclust below), and the method proposed in this paper, (referred to as K-H below). The biggest difference between the two methods is in the lacking ability of K-H to classify Reed warblers correctly. As for adult birds, the proposed method has better accuracy for Paddyfield warblers, and contrary to adult birds, also for Blyth’s reed warbler, although the differences are small.

	True species
	Reed warbler		Blyth’s reed w.		Paddyfield w.		Marsh warbler
Predicted species	mclust	K-H	mclust	K-H	mclust	K-H	mclust	K-H
Reed warbler	403	382	0	0	0	0	5	5
Blyth’s reed warbler	0	2	40	41	2	1	0	0
Paddyfield warbler	0	0	1	0	16	17	0	0
Marsh warbler	6	25	0	0	0	0	409	409

Table 8. Table 8: Quantiles of the posterior distribution for the parameters of the Acrocephalus model. The values in the ad -columns are the change in intercept for adult birds in the Regression parameter values -section, and separate values in the other sections. We have chosen to present the correlation between traits instead of the covariances. Notice the very slim posterior density intervals for Reed Warbler and Marsh Warbler , since we have lots of observations for these species, and the widely dispersed posterior distributions for Paddyfield Warbler, being a species with very few observations and a vague prior on the covariance matrices.

	Regression parameter values						Variance estimates						Correlation estimates
	Wing		Notch		Notch pos.		Wing		Notch		Notch pos.		Wing-N.		Wing-N.p.		N. - N.p.
Quantile	juv	ad	juv	ad	juv	ad	juv	ad	juv	ad	juv	ad	juv	ad	juv	ad	juv	ad
Reed warbler 2.5%	66.69	0.7	10.99	1.27	107.87	1.95	2.31	2.61	0.57	0.55	1.28	1.3	0.46	0.33	-0.09	0.00	0.46	0.29
50%	66.7	0.73	11.07	1.38	108	2.29	2.35	2.67	0.63	0.62	1.48	1.59	0.52	0.41	0.02	0.20	0.54	0.45
97.5%	66.72	0.76	11.13	1.52	108.13	2.56	2.38	2.74	0.71	0.72	1.66	2.35	0.57	0.49	0.13	0.37	0.61	0.60
Marsh warbler 2.5%	69.99	0.57	9.42	0.54	104.86	0.39	2.07	2.28	0.41	0.47	0.69	0.5	0.37	0.36	-0.13	-0.12	0.23	0.03
50%	70.05	0.69	9.45	0.61	104.96	0.66	2.19	2.51	0.44	0.53	0.8	0.76	0.40	0.43	-0.03	0.16	0.33	0.28
97.5%	70.1	0.81	9.48	0.68	105.04	1	2.32	2.76	0.47	0.6	0.92	1.16	0.44	0.49	0.06	0.41	0.42	0.49
Paddyfield warbler 2.5%	56.64	-0.74	12.24	0.4	114.8	0.05	1.28	2.36	0.49	0.76	0.77	0.68	-0.20	-0.08	-0.58	-0.11	-0.13	-0.03
50%	57.27	0.36	12.66	1.16	115.32	0.81	2.17	4.32	0.83	1.35	1.3	1.24	0.19	0.37	-0.25	0.35	0.27	0.42
97.5%	57.89	1.47	13.1	1.89	115.83	1.61	4	9.07	1.52	2.77	2.46	2.6	0.53	0.69	0.15	0.68	0.59	0.72
Blyth’s reed warbler 2.5%	61.84	-0.48	12.23	0.8	113.23	0.63	1.31	1.34	0.47	0.61	0.6	0.87	0.07	0.29	-0.18	-0.17	-0.09	0.09
50%	62.26	0.03	12.49	1.13	113.53	1.03	1.94	1.85	0.69	0.83	0.9	1.2	0.35	0.49	0.12	0.07	0.21	0.32
97.5%	62.69	0.56	12.75	1.46	113.85	1.43	3.04	2.64	1.08	1.17	1.37	1.71	0.58	0.64	0.40	0.30	0.47	0.51

Equations

$$Y_{ij}=(Y_{ij1},\dots ,Y_{ijq})\sim \text{N}\left(m_{ij},{\mathbf{\Sigma }}_{ij}\right)$$

$$m_{ij}=(m_{ij1},\dots ,m_{ijq})\text{and}{\mathbf{\Sigma }}_{ij}={\left({\Sigma }_{ijkl}\right)}_{k,l=1}^q$$

$$x_{ij}=\left(1,x_{ij1},\dots ,x_{ijp}\right)={\left(x_{ijm}\right)}_{m=0}^p$$

$${\mathbf{Y}}_i={\mathbf{X}}_i{\mathbf{B}}_i+{\mathbf{E}}_i$$

$${\mathbf{U}}_i=\text{vec}({\mathbf{Y}}_i)={\mathbf{Z}}_i{\beta }_i+\text{vec}({\mathbf{E}}_i)$$

$${\mathbf{Z}}_i={\mathbf{I}}_q\otimes {\mathbf{X}}_i=\left(\begin{array}{cccc}{\mathbf{X}}_i& 0& \cdots & 0\\ 0& {\mathbf{X}}_i& \ddots & ⋮\\ ⋮& \ddots & \ddots & 0\\ 0& \cdots & 0& {\mathbf{X}}_i\end{array}\right)$$

$${\mathbf{\Sigma }}_{ij}={\mathbf{\Sigma }}_i^{\alpha }\text{if }{x}_{ij}\in {\mathcal{X}}^{\alpha },$$

$$p({\beta }_i,{\mathbf{\Sigma }}_i^1,\dots ,{\mathbf{\Sigma }}_i^A)=p({\beta }_i)\prod _{\alpha =1}^{A}p({\mathbf{\Sigma }}_i^{\alpha })$$

$${\hat{\theta }}_i^{\text{(Bayes)}}=\frac{1}{R_i}\sum _{r=1}^{R_i}{\theta }_{ir}$$

$${\hat{\omega }}_i=\frac{1}{R_i}\sum _{r=1}^{R_i}f(Y;x,{\theta }_{ir})$$

$$Y_{ij}\in {\mathsf{S}}_{ij}={\mathsf{S}}_{ij1}\times \cdots \times {\mathsf{S}}_{ijq},$$

$${\mathsf{S}}_{ijk}=\{\begin{array}{ll}{Y}_{ijk},& k\notin {\mathsf{K}}_{ij},\\ \left(c_{ijk},d_{ijk}\right],& k\in {\mathsf{K}}_{ij},\end{array}$$

$${\mathsf{S}}_{ijk}=\{\begin{array}{ll}\left(-\infty ,\frac{1}{2}\right],& \text{if }{z}_{ijk}=0,\\ \left(z_{ijk}-\frac{1}{2},z_{ijk}+\frac{1}{2}\right],& \text{if }1\le z_{ijk}\le g_k-2,\\ \left(g_k-\frac{3}{2},\infty \right),& \text{if }{z}_{ijk}=g_k-1.\end{array}$$

$$Z_{ijk}=z({\mathsf{S}}_{ijk})=\{\begin{array}{ll}0,& \text{if }{\mathsf{S}}_{ijk}=\left(-\infty ,\frac{1}{2}\right],\\ \frac{c_{ijk}+d_{ijk}}{2},& \text{if Sijk is bounded},\\ g_k-1,& \text{if }{\mathsf{S}}_{ijk}=\left(g_k-\frac{3}{2},\infty \right],\end{array}$$

$${\mathcal{D}}^{\text{obs}}=\left\{\left(x_{ij},{\mathsf{S}}_{ij}\right);i=1,\dots ,N,j=1,\dots ,n_i\right\}.$$

$${\left\{{\theta }_{ir},Y_{ijkr},1\le j\le n_i,k\in {\mathsf{K}}_{ij};Y_{kr},k\in \mathsf{K}\right\}}_{r=1}^{R_i}$$

$$p({\theta }_i|{\mathcal{D}}_i^{\text{obs}})\prod _{j=1}^{n_i}f\left(y_{ij{\mathsf{K}}_{ij}}\mid x_{ij},{\mathsf{S}}_{ij},{\theta }_i\right)f\left(y_{\mathsf{K}}\mid x,Y_{{\mathsf{K}}^{\complement }};{\theta }_i\right)$$

$${\hat{\mathbf{\Sigma }}}_1^{\text{juv}}$$

$${\hat{\mathbf{\Sigma }}}_2^{\text{juv}}$$

$${\hat{\mathbf{\Sigma }}}_3^{\text{juv}}$$

$${\hat{\mathbf{\Sigma }}}_4^{\text{juv}}$$

$${\mathcal{D}}^{\text{new}}=(x,\mathsf{S})$$

$${\omega }_i={\iint }_{\mathsf{S}}f(Y;x,{\theta }_i)\prod _{k\in \mathsf{K}}\mathrm{d}{y}_{k}p({\theta }_i\mid {\mathcal{D}}_i^{\text{obs}})\mathrm{d}{\theta }_i,$$

$${\hat{p}}_i=\hat{\mathbb{P}}(I=i\mid {\mathcal{D}}^{\text{new}},{\mathcal{D}}^{\text{obs}})=\frac{{\pi }_i{\hat{\omega }}_i}{{\pi }_1{\hat{\omega }}_1+\dots +{\pi }_N{\hat{\omega }}_N},$$

$$\hat{\mathrm{I}}$$

$$=\underset{\mathcal{I}}{\mathrm{argmax}}\sum _{i=1}^{N}R(\mathcal{I},i)p_i$$

$${\Omega }_{\mathcal{I}}=\{(x,Y);\hat{\mathrm{I}}=\mathcal{I}\}$$

$$\Lambda =\bigcup _{|\mathcal{I}|>1}{\Omega }_{\mathcal{I}},$$

$$R(\mathcal{I},i)=\{\begin{array}{ll}0;& i\notin \mathcal{I}\\ 1/|\mathcal{I}|;& i\in \mathcal{I}\end{array}$$

$$\mathbb{E}\left[R(\mathcal{I},\mathrm{I})\mid {\mathcal{D}}^{\text{obs}},{\mathcal{D}}^{\text{new}}\right]=\frac{1}{|\mathcal{I}|}\sum _{i\in \mathcal{I}}{p}_i$$

$$\hat{\mathrm{I}}=\{(N)\}=\underset{i}{\mathrm{argmax}}{\pi }_i{\omega }_i$$

$$\begin{array}{cc}\hat{\mathbb{P}}\left(\hat{\mathrm{I}}\ne \mathrm{I}\mid x_1=0\right)=0.0251,& \hat{\mathbb{P}}\left(\hat{\mathrm{I}}\ne \mathrm{I}\mid x_1=1\right)=\mathrm{0.0264.}\end{array}$$

Full text

1 Introduction

1.1 Biological taxonomy

Biological taxonomy relies on correct classification of organisms to taxa, which in turn requires that reliable and detailed data is available to learn the decision rules. The use of traits to numerically classify taxa has a long history, cf. Sneath and Sokal (1973); Barnett et al. (1979); Payne and Preece (1980); Felsenstein (1983) for a summary of early work. These traits are either morphological or genetic, and measured on a scale that is either continuous, ordinal or categorical, or a mixture of these. In order to classify new organisms, some kind of reference material (or a training data set) is needed, to which the trait vector of the new subject is compared and matched. This requires that supervised learning or discriminant analysis has been applied to training data in order to generate a decision rule, so that the new organisms can be assigned to the taxon whose decision region includes this subject’s trait vector. Such a classifier is either phylogenetic or individual, depending on whether genetic relationships between species are taken into account or not when the decision regions are derived. For instance, in the context of bacterial classification, Wang et al. (2007) introduce an individual (naive Bayesian) classifier, whereas Gao et al. (2017) define a phylogenetic classifier such that when multiple taxa of the reference material match the trait vector of a new observation, information from a phylogenetic tree (Huelsenbeck and Ronquist, 2001) is used in order to facilitate classification.

Population ecology and ornithology in particular are important applications of biological taxonomy. To the authors’ knowledge, there has been no published interpretable and general approach to classifying ornithological data. The widely used work of Svensson (1992) in bird species identification makes rich use of numerical combinations (often refered to as “wing formulae”) of trait observations to classify birds to taxa, but offers no general solution. Modern approaches to bird species classification include the Merlin Bird ID software, by the Cornell Lab of Ornithology111See merlin.allaboutbirds.org.. Merlin Bird ID uses convolutional neural networks for image recognition (LeCun et al., 1989) and gives the user a ranked list of species to choose from. Due to the machine learning approach parameter interpretability is complicated though, and not accessible to the user.

1.2 Unsupervised learning

A number of unsupervised leaning methods methods have been proposed, as summarized for instance in Ripley (1996); Hastie et al. (2009); Liu et al. (2011). For continuous traits this includes linear (or Fisher) discriminant analysis (LDA, Fisher, 1936), quadratic discriminant analysis (QDA, Smith, 1946; Bensmail and Celeux, 1996; Fraley and Raftery, 2002), heteroscedastic discriminant analysis (Kumar and Andreou, 1998), regularized discriminant analysis (Friedman, 1989), mixture discriminant analysis (Hastie and Tibshirani, 1996), and probabilistic Fisher discriminant analysis (Bouveyron and Brunet, 2012). These are all soft discrimination methods (Wahba, 1998), where a joint statistical distribution of class membership (taxon) and trait is assumed, and the class with maximal posterior probability is selected. Bayes’ rule is used to compute the posterior probabilities in terms of prior probabilities of classes and the likelihoods of trait vectors conditional on classes. In particular, a naive Bayes classifier corresponds to a model for which the trait vectors are conditionally independent conditional on class (Domingos and Pazzani, 1997). Logistic regression and penalized logistic regression (Lin et al., 2000) can also be used for soft binary classification by specifying the posterior probabilities directly in terms of a regression model, with traits as covariates. In particular, this makes it easy to incorporate continuous as well as categorical traits. For any of the abovementioned soft classifiers, if the parameters of the likelihood are random as well, the resulting method is an instance of predictive classification or Bayesian discriminant analysis (Section 2.4 of Ripley, 1996; Aitchison and Dunsmore, 1980; Geisser, 1993; Hjort, 1986). This includes, for instance, the Bayesian version of QDA (Geisser, 1964; Fraley and Raftery, 2007). On the other hand, many machine learning methods, such as support vector machines (Cortes and Vapnik, 1995), distance-weighted discrimination (Marron et al., 2007), and boosting (Friedman et al., 2000), require no statistical model at all, since the decision regions of the classifier are defined in terms of an algorithm that does not rely on any model for how data was generated.

As mentioned above, it is often the case that taxonomic data is a mixture of various types of observations (continuous, integer valued, ordered categorical). For instance, there might be partial observations (due to e.g. rounding or grouping) or some components of the trait vector might be missing. Models for mixed type of data are important for a wide range of applications (De Leon and Chough, 2013). Here we will use the umbrella term obfuscated observations for partial and missing observations among the components of the trait vector. A general way to handle obfuscated data is to define a distance between pairs of trait vectors, and then employ a nearest neighbour method for classfication (Gower, 1971; Gower and Legendre, 1986; Kauffmann and Rousseeuw, 1990; D’Orazio, 2021). Another option for mixed data is to use a model-based soft classifier, with latent multivariate Gaussian variables, whose components are suiably obfuscated in order to generate the observed trait vector. These types of models have been used in Item Response Theory (IRT), originally proposed by Thurstone (1925) in the context of ecudational testing. The theory has since then been developed by many authors (Lord, 1952; Rasch, 1993; Vermunt, 2001; Lord and Novick, 2008). More recent Bayesian models of categorical data include Albert and Chib (1993); Fokoué and Titterington (2003); Chu and Ghahramani (2005); McParland et al. (2014, 2017). Another Bayesian soft classifier for ordinal data, a mixed membership model with a gamma distributed hidden variable, has been proposed by Virtanen and Girolami (2015).

It is often the case that two or more taxa fit the new observation almost equally well. It is possible then to include a reject option of the classifier that corresponds to “don’t know”. Binary classfication with reject options has been studied by Chow (1970) and Herbei and Wegkamp (2006) (see also Freund et al. (2004)), whereas Ripley (1996) adds a reject option for discrimination between an arbitrary number of classes. In taxonomy it is also important to provide some kind of cutoff when none of the candidate taxa fit the trait vector of the new organism well. This is achieved, for model-based soft classifiers, by declaring the newly observed trait vector as an outlier when its $p$-value is smaller than a cutoff. This $p$-value is based on a measure of oytlyingness of the observed new trait vector, using the mixture distribution of traits from training data of all classes as the analogue of a null distribution (Ripley, 1996).

1.3 The present article

The classification method presented in this paper is a generalization of the approaches suggested in Svensson (1992) and Malmhagen et al. (2013). We develop a method that makes it possible to conclude which trait measurements that characterize various taxa, as an alternative to neural networks and other black box methods. Our method can be viewed as an extension of Bayesian QDA that incorporates traits of mixed types, which are components of a hidden Gaussian vector that have been obfuscated in order to obtain the observed trait vector. In particular we draw on work by Albert and Chib (1993); Fokoué and Titterington (2003); McParland et al. (2014, 2017) and defined a Bayesian classifier whose parameters are estimated from training data by means of blockwise Gibbs sampling (Geman and Geman, 1984; Robert and Casella, 2013). The inclusion of covariates is mentioned as an extension in the discussion sections of McParland et al. (2014, 2017). Here we contribute to this extension, assuming that the location as well as the scale parameter of the latent Gaussian distribution depend on covariates, so that heteroscedasticity is achieved. We believe this covariate extension is important, since taxonomic trait measurments of population ecology data is usually informative conditional on other information. In other words, the same value of a trait measurement might indicate different taxa for a subject, depending on the particular subject’s additional information. In order to control for this additional information, we include it as set of covariates, and construct our classifier conditional on the covariate information for each subject.

A major novelty of our work is the use of set-valued outputs of the classifier. A set of size 0 corresponds to a scenario where no taxon fits the new organism (an outlier), a set of size 1 corresponds the case where one taxon fits the new organism better than all the others, whereas a set of size $N$, where $N$ is the number of classes, corresponds to a reject option. We also allow for decision regions with $2,\ldots,N-1$ species, corresponding to a partial reject option, where no single species is chosen, but some are excluded. In this way we generalize work by Chow (1970); Ripley (1996); Freund et al. (2004); Herbei and Wegkamp (2006), making use of a reward function with a set-valued input argument, so that an associated set-valued Bayesian classifier is defined. In order to handle outlying trait vectors, in contrast to Ripley (1996), a new organism’s trait vector is assigned one $p$-value for each taxon, rather than one single $p$-value for the whole training data set.

The primary scenario for which we envision use of our proposed method is as follows: A researcher has a desire to construct a procedure to classify subjects into a, usually small, number of taxa. The procedure should not require any advanced or expensive technology, but rather be viable using only visual observations or simple measurements from everyday tools, such as rulers, callipers and scales. Such a procedure has the potential to aid data collection for population ecologists, if high requirements on classification correctness are met. It is assumed that the researcher has access to training data where each subject has a know taxon, and the traits shared by the taxa are observed for the subjects, although not neccessarily perfectly. To achieve this goal, the researcher would apply our proposed methodology and, if feasible, derive a comprehensive list of which trait observations that signify the various taxa, and which are ambiguous. In short, the goal is to specify a decision rule that is not black box, which will be feasible when the number of traits is small.

The secondary scenario for using our methodology is that in which an agent desires to classify an organism to a taxon, but is unable to do so without some kind of aid. This could be due to inexperience of the considered taxa, or because the number of traits and the number of different outcomes of the traits are too large to extract a decision rule that is nimble. In such a situation, our method can provide tailored prediction of a taxon given the observed traits.

Our paper is organized as follows. In Section 2 we present the application with original ornithological data, which will be employed as we present and analyze our method of classification. The core models, multivariate and multiple regression with Gaussian errors for perfectly observed data, and the corresponding latent Gaussian models for obfuscated data, are presented in Section 3 in general but concise terms, and specified for our type of data. Section 4 presents soft (posterior probability) classification in conjunction with reward functions with set-valued input arguments. In Section 5 we propose model selection (both for traits and covariates) through cross-validation and apply it on another original data set. A comparison with mclust (Scrucca et al., 2016), the most popular R-package for clustering and classification using discriminant analysis, is made in Section 6. Section 7 contains the discussion. Most of the mathematical details are placed in the appendices, and we will refer to them accordingly.

As a final note, we stress that the presented approach in no way is limited to population ecology data and classification of subjects to taxa. It is applicable more generally for classification with partial reject options, of traits vectors whose components might be a mixture of different types.

2 Bird species classification

In order to showcase in detail the use of our method, we will present an example of bird species classification using an original data set A minor subset of this primary data set was used in Walinder et al. (1988) and Malmhagen et al. (2013), but the data set as a whole is unpublished. We will also make use of a second original, unpublished data set in Section 5, the details of which we present in that section. Both data sets were collected by the Falsterbo Bird Observatory.

The primary data set concerns four morphologically similar warblers in the Acrocephalus genus; Eurasian Reed Warbler (Acrocephalus scirpaceus), Marsh Warbler (Acrocephalus palustris), Blyth’s Reed Warbler (Acrocephalus dumetorum) and Paddyfield Warbler (Acrocephalus agricola). They constitute a typical example of species that experts with long experience can classify through visual observation, whereas those with less experience usually need measurements of particular traits to aid their classification. Our material contains measurements of these traits for birds that have been identified by experts.

Data was collected from wild birds captured and ringed in Falsterbo, Sweden, and from museum specimens at the museums of natural history in Stockholm, Copenhagen, and Tring. All birds have three traits of main concern, the wing length (measured as described in Svensson, 1992), the notch length of feather P2 (measured as described in Figure 1 of Malmhagen et al., 2013) and the relative position of the P2 notch to the rest of the wing, referred to as the notch position (taken as described in Svensson, 1992). These are all measurements of different parts of the wing. The only covariate is age, and it is binary, with levels juvenile and adult. Age is included as a covariate due to a suspected change in the trait distribution (both location and scale) between the birds’ juvenile and adult plumages.

Ideally, *wing length *and *notch length *are measured continuously, but in reality this is impossible. Instead, *wing length *is rounded to integer millimeters, as this is the typical unit that gives consistent measurements. By the same rationale *notch length *is rounded to half millimeters. Finally, *notch position *is by definition ordered categorical. It measures where in relation to the other wing feathers the notch is positioned and it is defined as the feather closest to the notch.

Overall we have 54 155 observed birds, and these constitute the referential data for this classification problem. These are distributed as presented in Table 1. The uneven distribution of observations across taxa will be commented on throughout the analysis of the data set.

3 Model formulation

We will first present the model for the case where all trait measurements are continuous and perfectly observed and thereafter for the case where the vector of trait measurements are obfuscated in different ways, including missing values.

3.1 Ideal case; no obfuscated trait measurements

Suppose we have $N$ different categories (or classes), contained in the set $\mathsf{N}=\{1,\ldots,N\}$, with known prior probabilities $\pi=(\pi_{1},\ldots,\pi_{N})$. With full data we measure $q$ traits and $p$ covariates of each subject. Let $Y_{ijk}$ be the measurement of trait $k$ for subject $j$ in category $i$, where $1\leq i\leq N$, $1\leq j\leq n_{i}$, $1\leq k\leq q$ and $n_{i}$ is the number of subjects in category $i$. We assume that

[TABLE]

are independent random vectors, each one having a multivariate normal distribution, with

[TABLE]

being the mean vector and the covariance matrix of subject $j$ of category $i$. The Gaussian assumption on $Y_{ij}$ is especially suitable when the inheritable components of the traits are known to be or can be assumed to be of polygenic nature and the environmental influence on the traits consists of many factors, each with a small effect. This is commonly the case for traits such as height or color of body parts (Lynch and Walsh, 1998; Lande et al., 2003). Let also

[TABLE]

be the covariate vector of subject $j$ of category $i$. Trait vectors and covariate vectors of category $i$ are rows in the matrices $\mathbf{Y}_{i}=\left(Y_{i1}^{\top},\ldots,Y_{in_{i}}^{\top}\right)^{\top}$ and $\mathbf{X}_{i}=\left(x_{i1}^{\top},\ldots,x_{in_{i}}^{\top}\right)^{\top}$ respectively, where $\top$ refers to matrix transposition. We now proceed by formulating a multivariate and multiple regression model

[TABLE]

for category $i$, where $\mathbf{B}_{i}=\left(B_{imk};m=0,\ldots,p;k=1,\ldots,q\right)$ is the regression parameter matrix, whose first row consists of intercepts for the $q$ traits, $m_{ij}$ is the $j^{\text{th}}$ row of $\mathbf{X}_{i}\mathbf{B}_{i}$, and $\mathbf{E}_{i}=\left(E_{i1}^{\top},\ldots,E_{in_{i}}^{\top}\right)^{\top}$ is an error term matrix with independent rows $E_{ij}\sim\text{N}(0,\mathbf{\mathbf{\Sigma}}_{ij})$.

For use in the construction of a joint prior, and later the derivation of the marginal posterior distributions of the parameters, the vectorized form of our regression model is needed. Denote by $\mbox{vec}(\cdot)$ the vectorization operation of a matrix, that appends the columns of the matrix from left to right on top of each other (Macedo and Oliveira, 2013). The inverse of the vec operation will also be used, and it is denoted by $\mbox{vec}^{-1}(\cdot)$. Then rewrite (2) as

[TABLE]

with $\beta_{i}=\text{vec}(\mathbf{B}_{i})$. Denoting an identity matrix of rank $q$ with $\mathbf{I}_{q}$ and using the matrix tensor (or Kronecker) product $\otimes$,

[TABLE]

is a block-diagonal matrix with $q$ blocks along the diagonal.

Now suppose we have $A$ covariance classes $\alpha=1,\ldots,A$ for category $i$ such that

[TABLE]

where $\mathcal{X}=\mathcal{X}^{1}\cup\ldots\cup\mathcal{X}^{A}$ is a disjoint decomposition of the predictor space $\mathcal{X}$. Assume a prior $\text{N}\left(\left(b_{i0k},\ldots,b_{ipk}\right)^{\top}=b_{ik},\mathbf{\Sigma}_{\mathbf{B}_{i}}\right)$ on each of the columns $k=1,\ldots,q$ of $\mathbf{B}_{i}$. This implies a prior $\text{N}\left(\left(b_{i1}^{\top},\ldots,b_{iq}^{\top}\right)^{\top}=\beta_{i0},\mathbf{I}_{q}\otimes\mathbf{\Sigma}_{\mathbf{B}_{i}}=\mathbf{\Sigma}_{\beta_{i}}\right)$ on $\beta_{i}$. Further, assuming prior independence and imposing an Inverse-Wishart distribution $\mathbf{\Sigma}_{i}^{\alpha}\sim IW(\nu_{0},\mathbf{V}_{0})$ (see for instance Gelman et al. (2013)) on the covariance matrices in (5) for $\alpha=1,\ldots,A$, we get the joint prior

[TABLE]

for the parameters of category $i$.

Write ${\cal D}_{i}=(\mathbf{X}_{i},\mathbf{Y}_{i})$ for the training data of category $i$ and let $(x,Y)$ be a new obervation that we want to classify. The predictive classifier in Section 4 will involve the density $\omega_{i}=f(Y;x|{\cal D}_{i})=\mathbb{E}[f(Y;x,\theta_{i}|{\cal D}_{i})]$ of this new observation, in case it belongs to category $i$, where $f(Y;x,\theta_{i})$ is the density function of the trait vector $Y$ conditional on the covariate vector $x$ and the parameter vector $\theta_{i}$. For a detailed derivation of the collection of model parameters $\theta_{i}=\left(\mathbf{B}_{i},\mathbf{\Sigma}^{1}_{i},\ldots,\mathbf{\Sigma}_{i}^{A}\right)$ and the posterior category weights $(\omega_{1},\ldots,\omega_{N})$ we refer to Appendix A. The Monte Carlo approximations of $\theta_{i}$ and $\omega_{i}$, $i\in\mathsf{N}$ are

[TABLE]

and

[TABLE]

respectively, where $R_{i}$ is the number of samples drawn from the posterior distribution of $\theta_{i}$, with $\theta_{ir}$ the parameter vector obtained by blockwise Gibbs sampling (Geman and Geman, 1984; Robert and Casella, 2013) in simulation run $r$.

3.2 General case; all types of obfuscation may occur

Overall our setup is the same as in Section 3.1, but now we suppose there is only partial information about the complete training data set ${\cal D}=\{(\mathbf{X}_{i},\mathbf{Y}_{i});\,i=1,\ldots,N\}$. Due to some obfuscation, which could be due to rounding, grouping, categorization or lost measurements of some traits, we typically have a mixture of different types of observations (De Leon and Chough, 2013) and only know that

[TABLE]

i.e. the complete trait vector $Y_{ij}$ for subject $j$ of category $i$ is contained in a hyperrectangle $\mathsf{S}_{ij}$, whose components are given by $\{\mathsf{S}_{ijk}\}_{k=1}^{q}$. These components are sets, ranging in possible size from singletons to infinite intervals of $\mathbb{R}$, and they are given by

[TABLE]

where $\mathsf{K}_{ij}=\left\{k;1\leq k\leq q;\,Y_{ijk}\text{ obfuscated}\right\}$. The obfuscations in (9)-(10) are of the same type as in IRT (Lord and Novick, 2008; Fox, 2010). The main difference is that the covariance matrix of the trait vector $Y_{ij}$ in (1) is arbitrary, whereas the covariance matrix of IRT models typically have a factor analysis structure (Spearman, 1904), so called mixture of factor analysis models for mixed data (MFA-MD; McParland et al., 2017).

More specifically, the obfuscations in (9)-(10) are of three main types. First, if a trait $Y_{ijk}$ is unobserved, written as $Y_{ijk}=\texttt{NA}$, the $k$:th component of $\mathsf{S}_{ij}$ is of infinite length; e.g. $c_{ijk}=-\infty$, $d_{ijk}=\infty$, and we let the interval be open. That is, the interval $\mathsf{S}_{ijk}$ equals $\mathbb{R}$. Secondly, a trait may be obfuscated in such a way that interval limits are observed. Rounding is a typical example of this; consider a measurement of a trait $y_{ijk}\in\mathbb{R}^{+}$ that has been rounded to $z_{ijk}\in 2\tau\cdot\mathbb{Z}^{+}$. We put $c_{ijk}=z_{ijk}-\tau$ and $d_{ijk}=z_{ijk}+\tau$, which constitute the limits of the interval around $z_{ijk}$. Generally, we assume rounding to the midpoint of an interval. We can always scale so that the interval is of unit length, which would be equivalent to $\tau=1/2$. Lastly we have the case when we observe an ordered categorical random variable $Z_{ijk}$. We assume there is an underlying normally distributed variable $Y_{ijk}$, and that each category $z_{ijk}$ corresponds to an interval of possible values of $y_{ijk}$. Count data can be treated as an instance of this type of obfuscation, for instance traits that represent the number of occurrences of something.

We will treat all types of obfuscations in the following unified way. Suppose trait $k$ of subject $j$ of category $i$ is imperfectly observed, i.e. $k\in\mathsf{K}_{ij}$. Let $g_{k}$ be the number of categories of this trait, which we number as $0,1,\ldots,g_{k}-1$. The observed category is $z_{ijk}\in\left\{0,1,\ldots,g_{k}-1\right\}$, where $g_{k}=2$ for binary data and $g_{k}=\infty$ for count data. The corresponding side of $\mathsf{S}_{ij}$ is

[TABLE]

Here, a useful trick would be to add auxiliary categories, that never were observed, to take the place of $z_{ijk}=0$ and $z_{ijk}=g_{k}-1$. That ensures all observed intervals are of unit length, although we may let intervals vary in length if there is reason to construct such a model. We also write

[TABLE]

for the center point of a finite or half-open and infinite $\mathsf{S}_{ijk}$, whereas $Z_{ijk}=z\left(\mathsf{S}_{ijk}\right)=Y_{ijk}$ when $Y_{ijk}=\mathsf{S}_{ijk}$ is perfectly observed. We will write the observed training data set as

[TABLE]

Finally, we remark on the importance (or lack thereof) of taking rounding into account. Consider rounding a Gaussian trait $Y_{ijk}$ to $z_{ijk}\in\mathbb{Z}$ for some $k$, $i=1,\ldots,N$ and $j=1,\ldots,n_{i}$. Suppose we have no covariates and that $Y_{ijk}\sim N(m_{ik},\sigma_{ik}^{2})$ for $j=1,\ldots,n_{i}$. An unbiased estimator of $m_{ik}$ is the average $\bar{Y}_{ik}$, whereas $\bar{Z}_{ik}$ is a biased estimator of $m_{ik}$. It is possible to quantify the size of the bias using $\sigma_{k}$ and the width of the rouding interval $\mathsf{w}=\eta_{k}\sigma_{k}$ (Tricker, 1984). In short, the larger $\mathsf{w}$ is relative to $\sigma_{k}$, the larger the bias is, as measured by $\eta_{k}$. Already when $\sigma_{k}=\mathsf{w}=1$, the bias is very small, and hence, unless an extremely precise mean estimate is needed, the bias is small compared to the uncertainty of the parameter estimate. Therefore, one might regard rounded values as true values, if the standard deviation of the trait that is rounded, is large enough.

Let $(x,\mathsf{S})$ refer to a new observation, for which the trait vector $Y\in\mathsf{S}$ is obfuscated for traits $k\in\mathsf{K}$. We refer to Appendices B and C for full details on the derivation of exact estimators of the model parameters $\theta_{i}$ and posterior weights $\omega_{i}=\mathbb{P}(\mathsf{S};x\mid\mathcal{D}_{i}^{\text{obs}})$ of category $i$, as well as Monte Carlo approximtions thereof. In brief, for each category $i$ we want to generate the set

[TABLE]

of $R_{i}$ blockwise Gibbs sampling iterates from the joint density

[TABLE]

of the parameters $\theta_{i}$, the imperfectly observed training data and the imperfectly observed new data point, where $\theta_{ir}=\left(\beta_{ir},\mathbf{\Sigma}_{ir}^{1},\ldots,\mathbf{\Sigma}_{ir}^{A}\right)$, $y_{ij\mathsf{K}_{ij}r}=\left(y_{ijkr};k\in\mathsf{K}_{ij}\right)$, and $y_{\mathsf{K}r}=\left(y_{kr};k\in\mathsf{K}\right)$ refer to the values of these quantities for Monte Carlo iteration $r$, whereas $Y_{\mathsf{K}^{\complement}}=\left(Y_{k};k\notin\mathsf{K}\right)$. In particular, $\theta_{ir}$ are drawn from the posterior distribution $p(\theta_{i}|{\cal D}^{\text{obs}}_{i})$ of $\theta_{i}$, conditional on observed data ${\cal D}^{\text{obs}}_{i}$ for category $i$. Whereas the Monte Carlo approximation of the Bayes estimator of $\theta_{i}$ is given by (7), also for obfuscated data, the Monte Carlo estimator of $\omega_{i}$ is slightly more complicated than (8), as described in Appendix B.2.

3.3 Example model fit

We now have the tools to fit a classification model to our Acrocephalus data described in Section 2. As we can tell from the end of Section 2, all traits are obfuscated in some way (many are missing), there is one covariate influencing the interpretation of the trait values and we have reason to believe there are different covariance classes.

Our trait vectors $Y_{ij}$ are attributed to species $i=1,\ldots,4$, where $i=1$ corresponds to Eurasian Reed Warbler, $i=2$ to Marsh Warbler, $i=3$ to Paddyfield Warbler and $i=4$ to Blyth’s Reed Warbler, whereas $j=1,\ldots,n_{i}$ denote individual birds within species. Each $Y_{ij}$ is of length $q=3$, where $Y_{ij1}$ is the wing length, $Y_{ij2}$ is the notch length and $Y_{ij3}$ is the notch position. All traits are obfuscated, but in different ways: $Y_{ij1}$ is continuous, but rounded to nearest integer on a millimeter scale; $Y_{ij2}$ is continuous, but rounded to nearest half millimeter; whereas $Y_{ij3}$ is ordered categorical.

We have one covariate age with values juveline or adult, which we code as 0 and 1 respectively, so that juvenile corresponds to the intercept. We denote this covariate by $x_{1}$, and it determines to which of the $A=2$ covariance classes each observation belongs. We denote the respective covariance matrices with $\mathbf{\Sigma}^{\texttt{juv}}$ and $\mathbf{\Sigma}^{\texttt{ad}}$.

The degrees of freedom hyperparameter of the covariance matrix $\mathbf{\Sigma}_{i}^{\alpha}$ of the prior is chosen as $\nu_{0}=10$, whereas the scale matrix $V_{0}$ has a diagonal of 15 and all other elements equal to 5. All covariance matrices have the same prior. Since there is $p=1$ covariate, the matrix used in the construction of the prior on the vectorized regression parameters, $\mathbf{\Sigma}_{\mathbf{B}_{i}}$, is a diagonal matrix with diagonal $(3,1)$ for all $i$. The mean parameter values $\mathbf{B}_{i0}=\mathbb{E}(\mathbf{B}_{i})$ of the prior on $\mathbf{B}_{i}$ are informative for each $i$ and based on the results in Malmhagen et al. (2013), as shown in Table 2.

Fitting the model in R (R Core Team, 2021), with all traits regarded as obfuscated, we get the Bayes estimates presented in Table 3. Highest posterior density intervals for each parameter are presented in Appendix D. Overall, the effect of age, our covariate, is to increase the trait values. However, the increase is different across traits and across species.

4 Classification

We now have interpretable results on the variation of traits, but ultimately we also want to use our knowledge to classify new birds to species. In this section we present classification for the model of Section 3.2. First, we define more generally the posterior category weights $\omega_{i}$ that were introduced in Section 3.1 and then look at canonical classification. Then we introduce set-valued classifiers, including the possibility of classifying empty sets in order to handle outliers. We also showcase the flexibility of the underlying method for classifying among a subset of categories and end with remarks on the choice of the classifier’s two tuning parameters $\rho$ and $\tau$.

Let

[TABLE]

denote a new observation with obfuscated traits $\mathsf{K}$. We define the posterior weight of category $i$ as

[TABLE]

where $f$ is the density function of the trait vector $Y=(y_{1},\ldots,y_{q})$ of the new observation, i.e. the multivariate Gaussian density function. As shown in Appendix B, the Markov chain in (11) can be used to find estimates $\hat{\omega}_{i}$ of these weights. We may then approximate the posterior probability $\hat{p}_{i}=\hat{\mathbb{P}}(I=i\mid{\cal D}^{\text{new}},{\cal D}^{\text{obs}})$ of ${\cal D}^{\text{new}}$ to be of category $i$ as

[TABLE]

where $I\in\mathsf{N}$ is the true but unknow category of the future observation, with prior distribution $\mathbb{P}(I=i)=\pi_{i}$.

Let $\mathcal{N}=\mathcal{P}(\mathsf{N})\setminus\emptyset$ denote the collection of all non-empty subsets of $\mathsf{N}$. Let $\hat{\mathrm{I}}\in\mathcal{N}$ be a classifier with $\lvert\hat{\mathrm{I}}\rvert\geq 1$. In order to define $\hat{\mathrm{I}}$ we introduce a reward function $\mathcal{N}\times\mathsf{N}\ni(\mathcal{I},i)\mapsto R(\mathcal{I},i)$ for all $\mathcal{I}\in\mathcal{N}$ and $i\in\mathsf{N}$. We interpret $R(\mathcal{I},i)$ as the reward of a classified set $\mathcal{I}$ of categories when the true category is $i$. Then put

[TABLE]

as the optimal predictive classifier or optimal Bayesian classifier, with the complete training data set ${\cal D}^{\text{obs}}$ and $p_{i}$ defined as in (15). Thus, $\hat{\mathrm{I}}$ is the set in $\mathcal{N}$ that maximizes the expected posterior reward. Each classifier $\hat{\mathrm{I}}=\hat{\mathrm{I}}({\cal D}^{\text{obs}},{\cal D}^{\text{new}})$, viewed as a function of a perfectly observed new data point ${\cal D}^{\text{new}}=(x,Y)$, partitions the test data space into decision regions

[TABLE]

for all $\mathcal{I}\in\mathcal{N}$. A similar definition of decision regions applies for an obfuscated new data point (13), where the trait part $\mathsf{S}$ is defined as rectangles of a grid, with $\Omega_{{\cal I}}$ the collection of $(x,\mathsf{S})$ for which $\hat{I}={\cal I}$. For any of these two definitions of decision regions, this gives rise to an indecisive region

[TABLE]

where we cannot distinguish one particular category with acceptable confidence, only eliminate some of the categories with low degree of belief.

Whenever the indecisive region $\Lambda$ is nonempty, the classifier allows for a partial reject option. An important special case when this occurs is Bayesian sequential analysis (DeGroot, 1970), for instance in the context of clinical trials (Carlin et al., 1998), before a final decision between the $N$ hypotheses (or categories) has been made. This corresponds to a partition of the test data space $\Omega$ into $N+1$ regions, so that either one of the $N$ hypotheses is chosen ($\hat{\mathrm{I}}=\{i\},i=1,\ldots,N$) or the decision is postponed to a later time point ($\hat{\mathrm{I}}=\{1,\ldots,N\}$).

There is considerable freedom in choosing the reward function $R$. For instance, Chow (1970) introduced a reward function for $N=2$ categories. This was generalized in Section 2.1 of Ripley (1996) to arbitrary $N\geq 2$, using $R(\{i\},\mathrm{I})=\mathbbm{1}(i\in\mathrm{I})$, $R(\{1,\ldots,N\},\mathrm{I})=c$ for some constant $0\leq c\leq 1$ and $R({\cal I},\mathrm{I})=0$ for all ${\cal I}$ with $2\leq|{\cal I}|\leq N-1$. In the next two subsections we will present two reward functions, the first of which corresponds to Ripley’s reward function when $|{\cal I}|=1$ and $|{\cal I}|=N$ (using $c=1/N$), but with a nonzero reward for $2\leq|{\cal I}|\leq N-1$.

4.1 Classification to one category

Let

[TABLE]

which has expected posterior reward

[TABLE]

and optimal classifier

[TABLE]

where $p_{(1)}<\ldots<p_{(N)}$ are the ordered posterior category probabilities. Notice that the indecisive region is empty ($\Lambda=\emptyset$), i.e. this reward function leads to a consistent, but potentially overzealous, classifier.

To estimate the probability of classifying wrongly using this approach, we simulated a large number of realisations $Y^{\star}$ from the predictive posterior distribution of the model in Section 3.3 for each species, under the assumption of a uniform prior distribution over species. Each $Y^{\star}$ was then attributed to a hyperrectangle in trait space in order to represent obfuscation as described in Section 2, resulting in half side lengths or $\tau$-values $(1/2,1/4,1/2)$ for trait $j=1,2,3$ respectively. We then computed which species an observation each hyperrectangle would predict, using the classifier in (17). This gives a numerical approximation under current obfuscations of the probability of observing a new bird and classifying it wrongly, when using the fitted model on the Acrocephalus data:

[TABLE]

Roughly 1 in 50 birds would be classified erraneously, when birds are distributed uniformly over the species under consideration. In Section 4.2 we will show how to reduce this error by allowing for classification to sets of species.

4.2 Classification with partial reject options

Choosing the reward function

[TABLE]

we get the expected posterior reward

[TABLE]

which is maximized by

[TABLE]

Thus we can tune the risk of classifying wrongly by picking $\rho\in[0,1]$ adequately, as it specifies an upper bound on the fraction of the largest posterior probability $p_{(N)}$ other posterior probabilities may attain and still be excluded. If we choose $\rho=0$, we get the classifier $\hat{\mathrm{I}}=\mathsf{N}$ which means $\mathbb{P}(\mathrm{I}\in\hat{\mathrm{I}})=1$ for all new observations, but that prediction method does not provide any information at all. The other extreme, choosing $\rho=1$, leads to $\hat{\mathrm{I}}=\{(N)\}$, and thus our classifier will be the same as (17). In Section 5.3 we present a way of choosing $\rho$ using cross-validation and a maximal accepted misclassification rate. In conclusion, our first classifier is a special case of the second.

Choosing $\rho=0.1$ yields the exclusion critera $\hat{p}_{i}<\hat{p}_{(N)}/10$. With this value of $\rho$, we find that the estimated probability of classifying wrongly using the Acrocephalus model rounded to four decimals are

[TABLE]

and that the probability of not singling out a particular species is

[TABLE]

This means we have reduced the probability of choosing the wrong species by $76.8\%$ for juvenile birds and $77.9\%$ for adult birds, at a price of not singling out a species with a probability of about 8% for any covariate value. Of the cases where $\lvert\hat{\mathrm{I}}\rvert>1$, only $0.0044\%$ will result in a classifier containing three species ($\lvert\hat{\mathrm{I}}\rvert=3$), meaning that we will be able to exclude at least half of the potential species for the vast majority of observations.

4.3 Classification with reject option

Following Section 2.1 of Ripley (1996), if none of the $N$ categories support test data ${\cal D}^{\text{new}}$ we would like to include $\emptyset$ as a possible output of the classifier $\hat{\mathrm{I}}$, so that $\hat{\mathrm{I}}\ \in\mathcal{P}(\mathsf{N})$. To this end, we denote the posterior weight of (14) as $\omega_{i}(x,\mathsf{S})$ in order to emphasize its dependence on the test data set (13). Then let

[TABLE]

where the outer integral is taken over all $\mathsf{S}^{\prime}$ such that $\omega_{i}(x,\mathsf{S}^{\prime})\leq\omega_{i}(x,\mathsf{S})$. We interpret $\bar{\omega}_{i}(x,\mathsf{S})$ as a $p$-value of test data $(x,\mathsf{S})$ for category $i$, i.e. the probability of observing an obfuscated trait vector $\mathsf{S}^{\prime}$ of category $i$ with covariate vector $x$, whose posterior weight $\omega_{i}(x,\mathsf{S}^{\prime})$ is at most as large as that of $(x,\mathsf{S})$. As such, it is a measure of the degree of outlyingness of ${\cal D}^{\text{new}}$. Our treatment differs from that of Ripley (1996) in that we define $N$ distinct $p$-values of test data, one for each category, whereas Ripley (1996) defines one single $p$-value for the mixture distribution of all categories. More specifically, given a value of $\rho$, we generalize the classifier (20) to

[TABLE]

where $\bar{\omega}_{i}=\bar{\omega}_{i}(x,\mathsf{S})$. Note that (20) is a special case of (22) with $\tau=0$.

Choosing $\tau=0.001$ results in

[TABLE]

The probability of not choosing a set containing the correct species is

[TABLE]

whereas the probability of not singling out a particular species ($|\hat{\mathrm{I}}|>1$) or getting an outlier ($|\hat{\mathrm{I}}|=0$) is

[TABLE]

These probabilities are very close to the ones in Section 4.2, meaning we can hedge the risk of classifying something that might be a new species (not belonging to $\{1,\ldots,N\}$) entirely at a low cost. The decision regions for these values on $\rho$ and $\tau$ are presented graphically in Appendix E, where we cover all classification scenarios with missing trait values as well.

It is also possible to include classification with empty outputs in the context of indifference zones (Bechhofer, 1954; Goldsman, 1986). Assume that the parameter space is divided into $N+1$ regions, the first $N$ of which correspond to each of the $N$ hypotheses (or categories), whereas the last region of the parameter space (the indifference zone) corresponds to scenarios where no hypothesis is adequate. Based on this, it is possible to divide the test data space into $N+1$ regions as well, depending on whether the posterior distribution of the parameter puts most of its probability mass in any of the first $N$ parameter regions ($\hat{\mathrm{I}}=\{i\},i=1,\ldots,N$) or in the indifference zone ($\hat{\mathrm{I}}=\emptyset$).

4.4 Subproblem accessibility

Having fitted the model to the whole set of species, one may use the fit for any subproblem, e.g. classifying a bird between two species when the others are, for some reason, ruled out. Taking species 1 and 2, i.e. Eurasian Reed Warbler and Marsh Warbler, we estimate the probability of classifying wrongly and the probability of ending up in the indecisive region $\Lambda$ analogously with Section 4.2. Using $\rho=0.1$ and $\tau=0$, we find that

[TABLE]

[TABLE]

for this particular subproblem.

4.5 Choosing $\tau$

A potential risk for misclassification is observing a subject of a category not even considered for classification. In order to mitigate this, we introduced $\tau$ as a cut-off value for the trait distributions. Indeed, the value of $\tau$ determines how large deviations in trait measurements we accept without suspecting that we actually observe a subject from an unconsidered category. Choosing $\tau=0$ allows us to classify any point in the whole trait space, i.e. we believe the model is perfect in the sense that no unconsidered categories will be observed.

5 Model selection using cross-validation

Model selection can be used to select covariates and/or traits from a larger set by trading goodness-of-fit against parsimony for each candidate model $m$. Bayesian model selection requires a prior distribution on $m$ (Smith and Spiegelhalter, 1980; Kass and Raftery, 1995; Green, 1995; Tadesse et al., 2005). This approach involves the computationally intractable likelihood ${\cal L}({\cal D}_{i}^{obs})$ for each category $i$ and candidate model $m$. Non-Bayesian approaches, on the other hand, typically involve the observed likelihood ${\cal L}(\hat{\theta}_{i};{\cal D}_{i}^{obs})$ for each category $i$ and model $m$, with $\hat{\theta}_{i}$ the ML-estimator of $\theta_{i}$, as well as another term that penalizes large models. This includes the AIC (Akaike, 1974, 1998), NIC (Murata, 1991), BIC (Schwarz et al., 1978; Fraley and Raftery, 1998), and integrated complete-data likelihood (ICL; Biernacki et al., 2000) criteria, as well as stepwise model selection procedures based on pairwise hypothesis testing (McLachlan and Rathnayake, 2014). Here we will use an intermediate approach that retains the Bayesian assumption of random $\theta_{i}$ for each model $m$ but then compare these models in terms of their predictive performance. To this end we will use $\kappa$-fold cross validation. This will be illustrated on another original, unpublished data set, on two subspecies of Common chiffchaff (Phylloscopus collybita), the collybita and abietinus subspecies. It contains measurements for birds classified to subspecies visually by an expert at Falsterbo Bird Observatory.

5.1 Cross-validation for our type of models

The idea of $\kappa$-fold cross-validation is well established, and used for a very wide range of model families, see e.g. Wood (2017, p. 256). Choosing $\kappa=n_{i}$ for category $i$ corresponds to the basic form of cross-validation (Cover, 1969; Stone, 1974). This procedure is however computationally expensive, because one has to fit as many models ($=\sum_{i=1}^{N}n_{i}$) as there are observations. Since the method under study is already computationally intensive, in particular under widespread obfuscation, large $q$ and large data sets, we recommend using $\kappa$-fold cross-validation with $\kappa$ a bit smaller, i.e. a non-exhaustive cross-validation. In an interesting paper Kohavi (1995) examines cross-validation in general when choosing between classifiers. Kohavi concludes that $\kappa<n_{i}$ is generally preferred when picking the best classifier using cross-validation.

To perform $\kappa$-fold cross-validation in general for our class of models, we begin by choosing $\kappa\in\mathbb{Z}^{+}$ independently of $i$. Then create fold $l$ for category $i$ by choosing uniformly at random a set $J_{il}\subset\{1,\ldots,n_{i}\}$ comprising $n_{i}/\kappa$ or $n_{i}/\kappa+1$ observations of ${\cal D}^{\text{obs}}_{i}$, the training data at hand for category $i$. Repeat this for all categories until we have left-out test data sets $J_{l}=\cup_{i=1}^{N}J_{il}$ for $l=1,\ldots,\kappa$. Then for each $l$ proceed to fit models on the observations ${\cal D}^{\text{obs}}_{(-l)}={\cal D}^{\text{obs}}\setminus\{(x_{ij},\mathsf{S}_{ij});j\in J_{il}\}_{i=1}^{N}$ that were not left out and then estimate the posterior category probabilities $\hat{p}_{1}({\cal D}^{\text{obs}}_{(-l)},(x_{ij},\mathsf{S}_{ij})),\ldots,\hat{p}_{N}({\cal D}^{\text{obs}}_{(-l)},(x_{ij},\mathsf{S}_{ij}))$ for each observation $(x_{ij},\mathsf{S}_{ij})\in J_{l}$ that was left out. Choosing a reward function $R$, the corresponding classifier $\hat{\rm{I}}$ may be applied to each set of posterior probabilities in order to generate predictions $\hat{\mathrm{I}}_{(-l)ij}$ of the category, or set of categories, of the left-out observations contained in $J_{l}$. As this is repeated for each $l$, all observations will have been left out at one point when the cross validation is finished.

To assess the predictive performance of the $M$ models under consideration, let $w_{i}>0$ be weights such that $\sum_{i=1}^{N}w_{i}=1$, and compute

[TABLE]

with a classifier (20) that corresponds to a prespecified value of $\rho$ and $\tau$, for $m=1,\ldots,M$, and an appropriately chosen reward function, such as (16) or (19). One could e.g. use the weights $w_{i}=n_{i}/\sum_{a=1}^{N}n_{a}$ or $w_{i}=1/N$, depending on whether it is more valuable to be able to predict a category with many observations or not. Based on (23), the best classifier is

[TABLE]

When having computed $R^{\texttt{cv}}_{m}$, for $m=1,\ldots,M$, one has the possibility to choose a more parsimonious model that performes almost as well as $\hat{m}$, under the usual tradeoff between simplicity and predictive performance.

5.2 Choosing traits and covariates

We will now examplify usage of cross validation with $\kappa=10$ folds for model selection. At hand we have a data set over two subspecies of Common chiffchaff (Phylloscopus collybita). Each subspecies has $q=9$ continuous traits measured. Although they are measured to half millimeters, we will regard them as perfect observations (${\cal D}^{\text{obs}}={\cal D}$), since the rounding bias is extremely small. Moreover we have two binary covariates (age and season), and we are interested in which combination of covariates and traits that best predicts the subspecies of a new bird. We will consider models with one covariate, both covariates and both covariates including an interaction term, meaning that the number of predictors $p$ ranges from 1 to 3.

For each covariate scenario, we fit a model $m$ with one trait, for each trait in turn, and then choose the one with the highest $R^{\texttt{cv}}$ value, based on the reward function (16). Keeping this trait, we then add another trait and repeat the procedure, and keep the two traits that yield the highest $R^{\texttt{cv}}$ value. The procedure is repeated until we have a model with 8 traits, and for each step, the value of $R^{\texttt{cv}}$ is stored. We used the weights $w_{i}=n_{i}/\sum_{k=1}^{N}n_{k}$ for all $i$, and put $\rho=1$ and $\tau=0$, i.e. used the classifier in (17) in the computation of $R^{\texttt{cv}}$.

The main reason for doing this forward selection-like procedure is to reduce the number of models that are fitted from about $(2^{9}-1)\cdot 4\cdot 10=20440$ (if every trait combination is tried) to $(9+8+\ldots+1)\cdot 4\cdot 10=1800$. Also, for each included trait, if the value of $R^{cv}$ does not increase, we may choose a more parsimonious model.

Figure 1 shows a plot of how $R^{\texttt{cv}}$ changes with the number of included traits for the various covariate combinations.

5.3 Choosing $\rho$

We will now present how the choice of $\rho$, for a classifier $\hat{I}$ based on a reward function (19), can be made more intuitive by letting the chosen value $\rho_{\delta}$ correspond to a highest acceptable misclassification rate $\delta$. To this end, we first introduce a reward function

[TABLE]

that outputs 1 when the true category belongs to the classified subset of categories, and otherwise 0. Interpreting $\tilde{R}(\mathcal{I},i)=1$ as a correct classification, we estimate the rate of correct classifications through cross-validation as

[TABLE]

in analogy with (23). In order to speed up computation, we first perform cross-validation to generate category weights $\omega_{1}(x_{ij},\mathsf{S}_{ij}),\ldots,\omega_{N}(x_{ij},\mathsf{S}_{ij})$ for each observaton $(x_{ij},\mathsf{S}_{ij})$ in ${\cal D}^{\text{obs}}$. We then compute the probabilities $p_{1}(x_{ij},\mathsf{S}_{ij}),\ldots,$ $p_{N}(x_{ij},\mathsf{S}_{ij})$ of all categories according to (15), given a prior distribution $(\pi_{1},\ldots,$ $\pi_{N})$ on the categories. The category probabilities can then be used to create classifiers $\hat{I}_{(-l)ij}$ for each observation $(x_{ij},\mathsf{S}_{ij})$ with $j\in J_{l}$, for a grid of $\rho$-values. These classifiers are then inserted into (24), alongside the correct category, and the misclassification rate $1-R_{\rho}^{cv}$ is computed. We may then pick

[TABLE]

To illustrate this, we performed $\kappa=31$-fold cross-validation for the acrocephalus data. Recall that the majority of observations in these data are partial, and thus classification is hard (see e.g. the decision regions in Appendix E). To compute the category weights, for each fold $J_{l}$ we used the Bayes estimates $\hat{\theta}_{(-l)i}$ as in (7), based on the observations $(x_{ij},\mathsf{S}_{ij})$ of category $i$ not left out ($j\notin J_{l}$), as parameter values for the trait distributions. Then we numerically integrated

[TABLE]

for each $(x_{ij},\mathsf{S}_{ij})$ with $j\in J_{l}$ and $l=1,\ldots,\kappa$. We used two different priors in the computation of $p_{1}(x_{ij},\mathsf{S}_{ij}),\ldots,p_{N}(x_{ij},\mathsf{S}_{ij})$; a uniform prior ($\pi_{i}=1/N$ for all $i$) and a prior proportional to the number of observations in each category ($\pi_{i}=n_{i}/\sum_{\iota=1}^{N}n_{\iota}$). As in Section 5.2, we chose $w_{i}=n_{i}/\sum_{k=1}^{N}n_{k}$ for all $i$. The classifiers were created according to (20) and for $\rho=0.01,0.02,\ldots,0.99$ we computed (25) using two different reward functions. First the reward function $\tilde{R}(\mathcal{I},i)$, defined in (24), and then a second reward function

[TABLE]

which is more restrictive, as we consider the classifier to be correct when it contains only the correct category. The second reward function (28) is not used for choosing $\rho$, but rather in order to illustrate how its misclassification rate differs from the one obtained from the other less restrictive reward function $\tilde{R}$. The misclassification rates as functions of $\rho$ are presented in Figure 2 for both priors. In Table 4 we give the misclassification rates for each category using the uniform prior, for $\rho=1$ and $\rho=0.3$.

The parameter $\rho$ can be used as a measure of how managable a classification problem is, given $\delta$. For any classification problem, we may first choose the misclassification rate upper limit $\delta$, and then find the highest value $\rho_{\delta}$ of $\rho$ for which the misclassification rate, based on reward function $\tilde{R}$, is below $\delta$. A more managable problem would then be one for which $\rho_{\delta}$ is higher. Classification problems where $\rho_{\delta}=1$ can be considered automizable, as there are sufficiently few dubious cases not to warrant any closer examination by an expert.

6 Comparison with discriminant analysis in the R package mclust

The mclust-package is the most popular package for multivariate clustering and classification (Scrucca et al., 2016) using R. It contains methods for clustering and classification of observations stemming from mixtures of Gaussian distributions, with cross-validation readily available to assess classification accuracy of the discriminant analysis models. We will use the discriminant analysis function MclustDA in order to compare its classification performance with our proposed method, using the Acrocephalus-data.

If one specifies prior distributions in the call to MclustDA, inference of model parameters is done through maximum aposteriori estimation (Fraley and Raftery, 2007). Then, in the next step these estimates are plugged into a Bayseian classifier with known parameters (Section 6.1 of Fraley and Raftery, 2002). In contrast, our classification scheme is fully Bayesian (14)-(15), with integration over posterior distributions of parameters. The priors on the model parameters differ somewhat, since we assume independent priors on the mean vector and covariance matrix of each class’ Gaussian distribution, whereas MclustDA does not (Fraley and Raftery, 2007). Moreover, we may specify priors for the mean vector for each class in our proposed method, while MclustDA does not permit this. Therefore, we chose the default option in MclustDA, which uses the mean of the traits from the data set. For our proposed method, we use the informative priors of Section 3.3. The prior on the species is uniform in the cross-validation scheme for both methods.

A further difference is that MclustDA do not permit the use of partial observations and covariates. Therefore, we reduce the Acrocephalus-data considerably to only consist of complete observations and partition it according to the binary covariate age. Consequently, two models are fitted with each method. These are then evaluated using 100-fold cross validation to estimate the prediction accuracy.

Results of the comparison of classification accuracy are presented in detail in Tables 5, 6 and 7. The take-away is that the methods preform similarily, except for juvenile birds, where our method to a larger extent than MclustDA classifies Reed warbler as Marsh warbler, and this is the source of its larger overall classification error probability among the juveniles.

Lastly, note also that the error probabilities with our method are higher than the estimated error probabilities (18) of Section 4.1. Although the estimates are computed in different ways, this nevertheless indicates the added value of including partial observations (as in Section 4.1) when available.

7 Discussion

Throughout this paper, we have defined and analysed a classification problem with two sets of observations for each subject; its trait vector $Y$ and its covariate vector $x$, where $x$ is informative about the interpretation of $Y$. Since the trait values are often subject to various types of obfuscation, we set up a unified Bayesian framework for these situations, using a latent multivariate Gaussian distribution, with parameters estimated through supervised learning and a blockwise Gibbs sampling algorithm. To formalize the classification, we introduced reward functions with a set-valued input argument and two tuning parameters $\rho\in[0,1]$ and $\tau\in[0,1)$. The choice of $\rho$ affects the size and location of the indecisive region $\Lambda$ of our discriminant rule with a partial reject option. This region is the part of observational space where our classifer does not have sufficient information to rule out all but one category, whereas $\tau$ puts a limit on how much we allow an observation to deviate from the bulk of the data and still allowing it to be classified by our decision rule. Finally, we present a method of covariate and/or trait selection, through cross-validation, in order to obtain classifers that are both parsimonious and efficient.

Overall, there are two main usages of the method presented in this paper. First, one may derive distinguishing characteristics (cf. Appendix E) of the categories considered, due to interpretability. Secondly, one may use a fitted model to classify new observations with statistical rigour. An example of the usefulness of the first case would be an ornithologist with a set of birds of known taxa, who doesn’t know what morphologically separates these taxa. Using this method, she may extract for which trait measurements there is a high probability of certain taxa and thereby create (and write down) an identification key. Further, if there are too many combinations of trait levels to memorize, the Bayesian procedure we have described may perform the classification in an automized way.

A number of extensions of our work are possible. The first one regards the distribution of the latent variable. Although a Gaussian distribution is frequently used within quantitative genetics, there are instances when other distributions are preferable. As trait distribution, due to the central limit theorem, the Gaussian is particularly suitable when the traits are assumed to be polygenic and possibly affected by environmental factors (Lynch and Walsh, 1998). A conceptually critical but often numerically negligable adjustment is to correct all latent Gaussian distributions by truncating them to positive values for some of the traits. The trait wing length in our first real world example has to be positive by definition, and hence we should adjust our Gaussian distribution accordingly. However, considering the values of the parameter estimates (see Table 3), it would essentially make no difference to impose such a restriction for this particular data set. In other cases, it could be more important. When multimodality and skewness is a concern, an extension to hidden variables with a normal mixture distribution is an option (Hastie and Tibshirani, 1996; McParland et al., 2014).

A second extension is to incluce nominal (unordered) traits. This can be achieved by incorporating elements of a multinomial probit model (Geweke et al., 1994; McParland et al., 2014). Although bird traits are typically ordered, it might be of interest to incorporate nominal traits when more diverse taxa are classified, or when genetic data in terms of unordered genotypes is part of the trait vector.

Third, it would be of interest to consider model selection methods other than cross-validation. For instance, a Monte Carlo approximation of the BIC criterion, the so called BIC-MCMC approach, performed well in some previous analyses (Frühwirth-Schnatter, 2011; McParland et al., 2017).

Fourth, the reliance on training data with known categories can potentially be relaxed, or at least partially so. Without a reference material, which serves as training data, the nature of our problem would correspond to unsupervised learning and to a large extent fall under general model-based clustering (Fraley and Raftery, 2002; McParland et al., 2014, 2017) of observations taken from a mixture of multivariate Gaussian distributions. To take a step towards unsupervised learning in our setting, one would need to allow for indecisive regions where classification is ambigous. Such a clustering algorithm with partial reject options would transfer ideas in Chow (1970); Ripley (1996); Herbei and Wegkamp (2006) and this paper of having incomplete decisions, to a framework of unsupervised learning. A challenge in this context is to incorporate the effect of covariates. Specifying the number of clusters in advance would make the problem more well behaved, but might oftentimes be impossible, since the fundamental problem to solve in practice would be to determine this number of clusters. Still, this would allow the method to be used in situations where it is not known how to classify observations at all, and thus make it possible to investigate how many clusters a given data set supports.

Fifth, modifying the method slightly in order to handle repeated measurements is a straightforward task within our multivariate Gaussian framework. The benefit with repeated measurements of the traits is a better understanding of the magnitude of the measurement error, when trait vectors of observations are replaced by averaged trait vectors, for all repeated measurements. One could then include the number of measurements into the classification method, with direct effects on the size of the indecisive region and the accuracy of the classifier.

Sixth, as mentioned in Section 3, we assume prior independence between the regression parameters of different categories. This allows the effect of a covariate to vary between the categories, as opposed to forcing apriori a covariate to have the same effect across categories. However, Appendix D lists the posterior means of the covariate effects from our example model of Section 3.3, and one may notice that the effect is similar for some traits across categories, and to some extent even across traits. This indicates a general effect of our covariate, and hence we could construct a model that emphasizes such a general effect, by introducing apriori dependencies between the regression parameters of different species.

Finally, it is of interest to apply and extend our discrminant analysis method with partial reject options in order to analyze data sets where some observations are known to belong to several clusters (Latouche et al., 2011). This requires an extended type of reward function, $R(\mathcal{I},\mathrm{I})$, where both the classification $\mathcal{I}$ and the true $\mathrm{I}$ are sets of categories.

Acknowledgements

The authors are grateful for the data on Acrocephalus warblers and Chiffchaffs provided by Falsterbo Bird Observatory, in particular Lennart Karlsson, Björn Malmhagen and Göran Walinder, and for helpful methodological feedback from Felix Wahl and Mathias Lindholm. We thank the Museum of Natural History in Stockholm, Copenhagen and Tring.

Appendix A Model formulation, complete data

Appendices A-E contain further mathematical details about the Bayesian model and classification procedure defined in Sections 3 and 4 respectively of the main article. In order to make the text self-contained, some details of the main article have been repeated. This Appendix A contains, in full detail, the derivation of estimators and the posterior category weights for a model using perfectly observed data. Posterior distributions for Bayesian multivariate linear regression with homoscedasticity assumption is readily available in Rossi et al. (2012, Sections 2.8 and 2.11), and we extend this to allow for heteroscedasticity.

Suppose we have $N$ different categories, contained in the set $\mathsf{N}=\{1,\ldots,N\}$, with prior probabilities $\pi=(\pi_{1},\ldots,\pi_{N})$. With full data we measure $q$ traits and $p$ covariates of each subject. Let $Y_{ijk}$ be the measurement of trait $k$ for subject $j$ in category $i$, where $1\leq i\leq N$, $1\leq j\leq n_{i}$, $1\leq k\leq q$ and $n_{i}$ is the number of subjects in category $i$. We assume that

[TABLE]

are independent random vectors having a multivariate normal distribution, with

[TABLE]

being the mean vector and the covariance matrix of subject $j$ of category $i$. Let also

[TABLE]

be the covariate vector of subject $j$ of category $i$. Trait vectors and covariate vectors of category $i$ are rows in the matrices $\mathbf{Y}_{i}=\left(Y_{i1}^{\top},\ldots,Y_{in_{i}}^{\top}\right)^{\top}$ and $\mathbf{X}_{i}=\left(x_{i1}^{\top},\ldots,x_{in_{i}}^{\top}\right)^{\top}$ respectively. We now proceed by formulating a multivariate and multiple regression model

[TABLE]

for category $i$, where $\mathbf{B}_{i}=\left(B_{imk};m=0,\ldots,p;k=1,\ldots,q\right)$ is the regression parameter matrix, whose first row consists of intercepts for the $q$ traits, $m_{ij}$ is the $j^{\text{th}}$ row of $\mathbf{X}_{i}\mathbf{B}_{i}$, and $\mathbf{E}_{i}=\left(E_{i1}^{\top},\ldots,E_{in_{i}}^{\top}\right)^{\top}$ is an error term matrix with independent rows $E_{ij}\sim\text{N}(0,\mathbf{\mathbf{\Sigma}}_{ij})$.

For use in the construction of a joint prior, and later the derivation of the marginal posterior distributions of the parameters, the vectorized form of our regression model is needed. Denote the operation of appending columns of a matrix by $\text{vec}(\cdot)$ (we will also use the inverse operation $\text{vec}^{-1}(\cdot)$ on column vectors) and rewrite (29) as

[TABLE]

with $\beta_{i}=\text{vec}(\mathbf{B}_{i})$. Denoting an identity matrix of rank $q$ with $\mathbf{I}_{q}$ and using the matrix tensor product $\otimes$,

[TABLE]

is a block-diagonal matrix with $q$ blocks along the diagonal.

Now suppose we have $A$ covariance classes $\alpha=1,\ldots,A$ for category $i$ such that

[TABLE]

where $\mathcal{X}=\mathcal{X}^{1}\cup\ldots\cup\mathcal{X}^{A}$ is a disjoint decomposition of the predictor space $\mathcal{X}$. Assuming a prior $\text{N}\left(\left(b_{i0k},\ldots,b_{ipk}\right)^{\top}=b_{ik},\mathbf{\Sigma}_{\mathbf{B}_{i}}\right)$ on each of the columns $k=1,\ldots,q$ of $\mathbf{B}_{i}$, we obtain a prior $\text{N}\left(\left(b_{i1}^{\top},\ldots,b_{iq}^{\top}\right)^{\top}=\beta_{i0},\mathbf{I}_{q}\otimes\mathbf{\Sigma}_{\mathbf{B}_{i}}=\mathbf{\Sigma}_{\beta_{i}}\right)$ on $\beta_{i}$. Further, assuming prior independence of $\beta_{i},\mathbf{\Sigma}_{i}^{1},\ldots,\mathbf{\Sigma}_{i}^{A}$ and imposing an Inverse-Wishart distribution $\mathbf{\Sigma}_{i}^{\alpha}\sim IW(\nu_{0},\mathbf{V}_{0})$ on the covariance matrices in (32) for $\alpha=1,\ldots,A$, we get the joint prior

[TABLE]

for the parameters of category $i$.

A.1 Estimation

Let $\theta_{i}=\left(\mathbf{B}_{i},\mathbf{\Sigma}^{1}_{i},\ldots,\mathbf{\Sigma}_{i}^{A}\right)$ represent all parameters of category $i$. In the following, we assume that $\theta_{1},\ldots,\theta_{N}$ are independent random vectors with probability densities $p(\theta_{1}),\ldots,p(\theta_{N})$ defined in (33). Introducing dependencies is of course possible, and may be important for specific problems. This is briefly mentioned in Section 7 of the main paper. From Bayes’ Theorem we get an aposteriori density

[TABLE]

of $\theta_{i}$ given the complete training data set ${\cal D}_{i}=\{(x_{ij},Y_{ij});\,j=1,\ldots,n_{i}\}=\{\mathbf{X}_{i},\mathbf{Y}_{i}\}$ for category $i$. The function ${\cal L}\left(\theta_{i};{\cal D}_{i}\right)=p\left(\mathbf{Y}_{i}\mid\mathbf{X}_{i},\theta_{i}\right)$ is the likelihood, whereas ${\cal L}({\cal D}_{i})=p({\cal D}_{i})$ is the marginal likelihood of category $i$. In the last step we removed the normalizing factor ${\cal L}({\cal D}_{i})^{-1}$, since it does not depend on $\theta_{i}$. The Maximum Aposteriori (MAP)-estimator of $\theta_{i}$ is

[TABLE]

whereas the Bayes’ estimator of $\theta_{i}$ is

[TABLE]

Finally, given a new observation ${\cal D}^{\text{new}}=(x,Y)$, define the posterior probability of the new observation belonging to category $i$ as

[TABLE]

where

[TABLE]

are the posterior category weights given ${\cal D}^{\text{new}}$ for all categories, before the prior probabilities $\pi_{i}$ have been taken into account.

A.2 Monte Carlo Approximations

It is usually difficult to evaluate the normalizing constants ${\cal L}({\cal D}_{i})^{-1}$ for high-dimensional data sets, and hence also $\theta^{(\text{Bayes})}_{i}$ and $\omega_{i}$. However, it is possible to estimate $\theta^{(\text{Bayes})}_{i}$ and $\omega_{i}$ by Monte Carlo simulation, with

[TABLE]

and

[TABLE]

respectively, if $\theta_{i1},\ldots,\theta_{iR_{i}}$ are $R_{i}$ replicates drawn from the posterior distribution $p(\theta_{i}\mid{\cal D}_{i})$, with $\theta_{ir}=\left(\beta_{ir},\mathbf{\Sigma}_{ir}^{1},\ldots,\mathbf{\Sigma}_{ir}^{A}\right)$.

We will generate $\theta_{i1},\ldots,\theta_{iR_{i}}$ by blockwise Gibbs sampling, and for this we need the conditional posterior distributions of $\beta_{i}$ and $\mathbf{\Sigma}_{i}^{\alpha}$ for $\alpha=1,\ldots,A$. To derive those, we need some additional notation. Let $\mathbf{Z}_{i}^{\alpha}$, $\mathbf{X}_{i}^{\alpha}$, $\mathbf{Y}_{i}^{\alpha}$ and $\mathbf{U}_{i}^{\alpha}$ denote the submatrices of $\mathbf{Z}_{i}$, $\mathbf{X}_{i}$, $\mathbf{Y}_{i}$ and $\mathbf{U}_{i}$ corresponding to covariance class $\alpha$, and let $\mathbf{I}_{n}$ be the identity matrix of order $n$. Recall also that $\mathbf{B}_{i}=\text{vec}^{-1}(\beta_{i})$, meaning that we know $\mathbf{B}_{i}$ from $\beta_{i}$, and vice versa. For simplicity of notation we omit index $i$ in the following proposition:

Proposition 1**.**

Denote the parameter vector of a Bayesian multivariate multiple regression model with $A$ covariance classes by $\theta=\left(\beta,\mathbf{\Sigma}^{1},\ldots,\mathbf{\Sigma}^{A}\right)$, where $\beta$ is the regression parameter vector and $\mathbf{\Sigma}^{1},\ldots,\mathbf{\Sigma}^{A}$ are the $A$ covariance matrices. Let the prior of $\theta$ be $p(\theta)=p(\beta)\prod_{\alpha=1}^{A}p(\mathbf{\Sigma}^{\alpha})$, where $\beta\sim\text{N}(\beta_{0},\mathbf{\Sigma}_{\beta})$ and $\mathbf{\Sigma}^{\alpha}\sim IW(\nu_{0},\mathbf{V}_{0})$ for $\alpha=1,\ldots,A$. Then the posterior distribution of $\beta\mid\mathbf{U},\mathbf{Z},\mathbf{\Sigma}^{1},\ldots,\mathbf{\Sigma}^{\alpha}$ is $\text{N}(\tilde{\beta},\tilde{\mathbf{\Sigma}})$, where

[TABLE]

and

[TABLE]

Proof.

By applying Bayes’ theorem

[TABLE]

where $n^{\alpha}$ is the number of observations in covariance class $\alpha$,

[TABLE]

where in the second step of the last equation we used Lemma 1 below, and

[TABLE]

Consequently,

[TABLE]

where

[TABLE]

and

[TABLE]

Notice that $\tilde{\beta}$ is a multivariate version of a weighted average of the prior vector $\beta_{0}$ and the least squares estimates of $\beta$ obtained from the $A$ covariance classes. ∎

The following lemma was used in the proof of Proposition 1. It admits a considerable gain in computational speed, when calculating the posterior covariance matrix $\tilde{\mathbf{\Sigma}}$.

Lemma 1**.**

Let $\mathbf{Z}$ be a block-diagonal $nq\times(p+1)q$ matrix, where there are $q$ blocks $\mathbf{X}$, which are $n\times(p+1)$-matrices, along the diagonal. Let $\mathbf{\Sigma}$ be a symmetric, positive definite $q\times q$-matrix. Then it holds that

[TABLE]

Proof.

We prove the lemma by iterated use of the mixed-product property of the tensor product. Since $\mathbf{Z}=\left(\mathbf{I}_{q}\otimes\mathbf{X}\right)$, the left hand side becomes

[TABLE]

and the right hand side becomes

[TABLE]

which proves the lemma already in the third equalities. ∎

Using the vector form, we may express the conditional posterior of the regression parameters

[TABLE]

where

[TABLE]

and

[TABLE]

Meanwhile, the conditional posteriors of the covariance matrices are

[TABLE]

where $n_{i}^{\alpha}$ denotes the number of observations in the covariance class $\alpha$ for category $i$ and

[TABLE]

Having computed $\hat{\omega}_{1},\ldots,\hat{\omega}_{N}$, for ${\cal D}^{\text{new}}$, we may compute the Monte Carlo-estimated aposteriori probability of ${\cal D}^{\text{new}}$ being in category $i$ as

[TABLE]

where ${\cal D}={\cal D}_{1}\cup\ldots\cup{\cal D}_{N}$ is the complete training data set.

Appendix B Model Formulation, obfuscated data

Overall our setup is the same as in Appendix A, but we now suppose there is only partial information about the complete training data set ${\cal D}$. Due to some obfuscation, which could be due to rounding, grouping, categorization or lost measurements of some traits, we only know that

[TABLE]

i.e. the complete trait vector $Y_{ij}$ for subject $j$ of category $i$ is contained in a hyperrectangle $\mathsf{S}_{ij}$, whose components are given by $\{\mathsf{S}_{ijk}\}_{k=1}^{q}$. The components are sets, ranging in possible size from singletons to infinite intervals of $\mathbb{R}$, and are given by

[TABLE]

where $\mathsf{K}_{ij}=\left\{k;1\leq k\leq q;\,Y_{ijk}\text{ obfuscated}\right\}$. As described in the main article, without loss of generality we assume that $z_{ijk}$, the mid point of $(c_{ijk},d_{ijk}]$ for finite sets, is integer-valued.

We can treat all types of obfuscations uniformly in the following way. Suppose trait $k$ of subject $j$ of category $i$ is imperfectly observed, i.e. $k\in\mathsf{K}_{ij}$. Let $g_{k}$ be the number of categories of this trait, which we number as $0,1,\ldots,g_{k}-1$. The observed category is $z_{ijk}\in\left\{0,1,\ldots,g_{k}-1\right\}$, where $g_{k}=2$ for binary data and $g_{k}=\infty$ for count data. The corresponding side of $\mathsf{S}_{ij}$ is

[TABLE]

We also write

[TABLE]

for the center point of a finite or half-open, infinite $\mathsf{S}_{ijk}$, whereas $z\left(\mathsf{S}_{ijk}\right)=Y_{ijk}$ when $Y_{ijk}=\mathsf{S}_{ijk}$ is perfectly observed. We may write the actually observed training data set as

[TABLE]

B.1 Estimation

Using ${\cal D}^{\text{obs}}_{i}=\{(x_{ij},\mathsf{S}_{ij});j=1,\ldots,n_{i}\}$, the posterior distribution of $\theta_{i}$ becomes

[TABLE]

where the normalizing factor

[TABLE]

is the marginal likelihood of category $i$, and

[TABLE]

Thus, with perfect observations ($\mathsf{K}_{ij}=\emptyset$), we evaluate the density of the trait vector $f$ at the observed point $Y_{ij}$ and the model is exactly as specified in Appendix A. Otherwise, we construct a $\left|\mathsf{K}_{ij}\right|$-dimensional integral over $f$ and the contribution to the likelihood is this integral, with the function evaluated exactly at the remaining perfectly observed traits, if any exist. In particular, if all traits are imperfectly observed, the integral is $q$-dimensional. We may approximate the integral in (38) by

[TABLE]

whenever all $\left|\mathsf{S}_{ijk}\right|<\infty$ for $k\in\mathsf{K}_{ij}$, which is the case when employing the trick with auxiliary categories, mentioned in Section 3.2 of the main article.

Alternatively, since $p(\mathsf{S}_{ij},x_{ij};\theta_{i})$ potentially contains integrals of a multivariate Gaussian density function and there in general is a lack of a CDF on closed form for this distribution, the integrals in (38) need to be solved numerically. However, in the case of $|\mathsf{K}_{ij}|=1$, with $\mathsf{K}_{ij}=\{k\}$ and $\mathsf{S}_{ijk}=\left(c_{ijk},d_{ijk}\right]$, the integral is univariate and thus222The notation with subscript $(-k)$ means dropping element $k$ from a vector; dropping row $k$ from a matrix when not being the last index of a matrix; and dropping column $k$ when being the last index.

[TABLE]

where $m_{ijk}(y_{ij(-k)})=m_{ijk}+\mathbf{\Sigma}_{ijk(-k)}\mathbf{\Sigma}_{ij(-k)(-k)}^{-1}(y_{ij(-k)}-m_{ij(-k)})$ is the conditional expectation of $Y_{ijk}$ given that $Y_{ij(-k)}=(Y_{ijk^{\prime}};\,k^{\prime}\neq k)=y_{ij(-k)}$, $\sigma_{ijk}=\sqrt{\mathbf{\Sigma}_{ijkk}-\mathbf{\Sigma}_{ijk(-k)}\mathbf{\Sigma}_{ij(-k)(-k)}^{-1}\mathbf{\Sigma}_{ij(-k)k}}$ is the conditional standard deviation of $Y_{ijk}$ given any value of $Y_{ij(-k)}$, and $\Phi$ is the CDF of the univariate Gaussian distribution with mean 0 and standard deviation 1.

Using ${\cal D}^{\text{obs}}_{i}$, we find from (B.1) that the estimators $\theta^{(\text{MAP})}_{i}$ and $\theta^{(\text{Bayes})}_{i}$ are

[TABLE]

and

[TABLE]

respectively. Furthermore, redefining ${\cal D}^{\text{new}}:=(x,\mathsf{S})$ for a new observation, where $\mathsf{S}=\mathsf{S}_{1}\times\ldots\times\mathsf{S}_{q}$, and denoting the corresponding set of imperfectly observed traits by $\mathsf{K}$, leads to the posterior category weights

[TABLE]

of this observation.

B.2 Monte Carlo Approximations

The integral over $\mathsf{S}$ in (41) is, as mentioned in conjunction with (39), potentially impossible to compute analytically, but could also be well behaved. We can in theory approximate $\theta^{(\text{Bayes})}_{i}$ in (B.1) as in (35) by sampling $\theta_{i}$ from $p(\theta_{i}\mid{\cal D}^{\text{obs}})$ a total of $R_{i}$ times. However, this entails a large number of numerical evaluations of integrals, see (B.1)-(38). Similarly, we may estimate $\omega_{i}$ for $1\leq i\leq N$ in (41) through

[TABLE]

which in addition to previously presented integrals, involves computation of an integral over $\mathsf{S}$. As an alternative way of computing (B.1) and (42), we also present an approach where complete data is sampled, based on the obfuscated data, as one step of the Monte Carlo algorithm, whereas the parameters are sampled as another step of the same algorithm. This allows us to estimate all $\theta^{(\text{Bayes})}_{i}$ and $\omega_{i}$ under widespread obfuscation, given that we are able to simulate $\mathbf{Y}_{i}$, $i=1,\ldots,N$. Overall, we want to generate

[TABLE]

from

[TABLE]

where $\theta_{ir}=\left(\beta_{ir},\mathbf{\Sigma}_{ir}^{1},\ldots,\mathbf{\Sigma}_{ir}^{A}\right)$, $y_{ij\mathsf{K}_{ij}r}=\left(y_{ijkr};k\in\mathsf{K}_{ij}\right)$, $y_{\mathsf{K}r}=\left(y_{kr};k\in\mathsf{K}\right)$ and $Y_{\mathsf{K}^{\complement}}=\left(Y_{k};k\notin\mathsf{K}\right)$. Note that we do not condition on $\mathsf{S}$ in the conditional density of the unobserved traits for the new observation that we want to classify, as this would introduce a bias in the Monte Carlo estimate of $\omega_{i}$ below.

The details of the specific Gibbs sampling approach we use are presented in Appendix C. Having generated a sample $\theta_{i1},\ldots,\theta_{iR_{i}}$, we may compute the estimated category weights of ${\cal D}^{\text{new}}$ as

[TABLE]

where $Y_{\mathsf{K}^{\complement}}$ is as above, and $Y_{\mathsf{K}r}=\left(Y_{kr};k\in\mathsf{K}\right)$. For every $\theta_{ir}$, one could generate many $Y_{\mathsf{K}}$ and replace the indicator with an average of the indicators for each sampled $Y_{\mathsf{K}}$.

A potentially more efficient method would be to define $\{y_{t}\}_{t=1}^{T}$, where $y_{t}=(y_{t1},\ldots,y_{tq})$ with $y_{t\mathsf{K}^{\complement}}=Y_{\mathsf{K}^{\complement}}$ and $y_{t\mathsf{K}}\in\bigtimes_{k\in\mathsf{K}}\mathsf{S}_{k}$, in such a way that $\{y_{tk},k\in\mathsf{K}\}$ is a grid approximation of $\bigtimes_{k\in\mathsf{K}}\mathsf{S}_{k}$. Then we can estimate $\omega_{i}$ through

[TABLE]

If we use the trick with auxiliary categories described in Section 3.2 of the main article, we can we can choose $y_{t}$ uniformly at random on $\bigtimes_{k\in\mathsf{K}}\mathsf{S}_{k}$, as long as we do not have any missing observations, since those are represented with an infinite interval. Thus, (46) is potentially more effcient than (45), but comes at a cost of generality, since (45) is applicable to any new observation.

Finally, the Monte Carlo-estimated aposteriori probability of ${\cal D}^{\text{new}}=(x,\mathsf{S})$ being of category $i$ is

[TABLE]

and we may apply (41) with replacement of $\omega_{i}$ by (45) for prediction. If (46) is inserted into (47) we notice that $\prod_{k\in\mathsf{K}}\lvert\mathsf{S}_{k}\rvert$ and $T$ cancel out, and in case $R_{i}=R$ for $i=1,\ldots,N$, also $R$ cancels out.

Appendix C Gibbs sampling details

The focus of this appendix is Procedure 1, in which we describe in detail how to generate a sample of size $R_{i}$ from the posterior distribution of the parameter vector $\theta_{i}$, using blockwise Gibbs sampling. It describes the general case, i.e. when we have obfuscated trait measurements in ${\cal D}^{\text{obs}}$. For the case with perfectly observed trait measurements, we skip the sampling of $\mathbf{Y}_{i}$ and use the observed values instead, otherwise the procedure is the same. Applying the procedure to data from each category $i$ will yield all the samples we need from the posterior distribution in order to perform classification.

In Procedure 1, $TN(\mu,\mathbf{\Sigma},\mathsf{S})$ refers to the truncated Gaussian distribution, where $\mu$ is the mean vector, $\mathbf{\Sigma}$ is the covariance matrix and $\mathsf{S}$ is a hyper rectangle specifying the truncation limits. Simulating from this distribution can be done exactly using rejection sampling, or approximately using an inner Gibbs algorithm. Depending on application, either approach can be preferred, as the tradeoff is exact sampling versus efficiency. Also, more advanced algorithms such as Importance Sampling-techniques can be used in this step.

Appendix D Highest posterior density intervals of parameters

Appendix E Visualized decision regions

All of these visualizations are done using the same model fit and the same generated new observations from the posterior predictive distribution as in Section 4 of the main text. As a reminder, we used the values $\rho=0.1$ and $\tau=0.001$ for the tuning parameters of the classifier.