An information criterion for auxiliary variable selection in incomplete data analysis

TL;DR

This paper introduces an information criterion designed to select relevant auxiliary variables in incomplete data analysis, improving primary variable estimation by effectively leveraging related auxiliary data.

Contribution

It proposes a novel information criterion that is asymptotically unbiased for model selection in incomplete data settings, linking it to leave-one-out cross validation.

Findings

The criterion effectively identifies relevant auxiliary variables in simulations.

It improves estimation accuracy for primary variables when auxiliary variables are relevant.

Performance demonstrated on real data example.

Abstract

Statistical inference is considered for variables of interest, called primary variables, when auxiliary variables are observed along with the primary variables. We consider the setting of incomplete data analysis, where some primary variables are not observed. Utilizing a parametric model of joint distribution of primary and auxiliary variables, it is possible to improve the estimation of parametric model for the primary variables when the auxiliary variables are closely related to the primary variables. However, the estimation accuracy reduces when the auxiliary variables are irrelevant to the primary variables. For selecting useful auxiliary variables, we formulate the problem as model selection, and propose an information criterion for predicting primary variables by leveraging auxiliary variables. The proposed information criterion is an asymptotically unbiased estimator of the Kullback-Leibler divergence for complete data of primary variables under some reasonable conditions. We also clarify an asymptotic equivalence between the proposed information criterion and a variant of leave-one-out cross validation. Performance of our method is demonstrated via a simulation study and a real data example.

Figures

Tables

Table 1. Table 1: Random variables in incomplete data analysis with auxiliary variables. B = ( Y , A ) 𝐵 𝑌 𝐴 B=(Y,A) is used for estimation of unknown parameters, and X = ( Y , Z ) 𝑋 𝑌 𝑍 X=(Y,Z) is used for evaluation of candidate models.

	Observed	Latent	Complete
Primary	$Y$	$Z$	$X$
Auxiliary	$A$	–	–
All	$B$	–	$C$

Table 2. Table 2: Comparisons between θ ^ b subscript ^ 𝜃 𝑏 \hat{\theta}_{b} and θ ^ y subscript ^ 𝜃 𝑦 \hat{\theta}_{y} for predicting X 𝑋 X , and that for Y 𝑌 Y .

	$p_x\left({\hat{\theta }}_b\right)$ vs. $p_x\left({\hat{\theta }}_y\right)$	$p_y\left({\hat{\theta }}_b\right)$ vs. $p_y\left({\hat{\theta }}_y\right)$
	${\mathrm{AIC}}_{x;b}-{\mathrm{AIC}}_{x;y}$	${\mathrm{AIC}}_{y;b}-{\mathrm{AIC}}_{y;y}$
Case 1	-2.67	-0.96
Case 2	9.86	10.37

Table 3. Table 3: Expected AIC difference is compared with the risk difference. The values are computed from T = 10 4 𝑇 superscript 10 4 T=10^{4} runs of simulation with their standard errors in parentheses.

$n$	100	200	500	1000	2000	5000
$E\left[{\mathrm{AIC}}_{x;b}-{\mathrm{AIC}}_{x;y}\right]$	-3.559	-3.263	-3.221	-3.197	-3.195	-3.180
	(0.074)	(0.021)	(0.015)	(0.013)	(0.013)	(0.012)
$2n\left\{{ℛ}_x\left({\hat{\theta }}_b\right)-{ℛ}_x\left({\hat{\theta }}_y\right)\right\}$	-3.603	-3.333	-3.275	-3.208	-3.182	-3.232
	(0.071)	(0.054)	(0.050)	(0.050)	(0.050)	(0.050)

Table 4. Table 4: Useful auxiliary variable (Case 1): selection frequencies of θ ^ b subscript ^ 𝜃 𝑏 \hat{\theta}_{b} and θ ^ y subscript ^ 𝜃 𝑦 \hat{\theta}_{y} .

$n$	100	200	500	1000	2000	5000
${\hat{\theta }}_b$	9230	9475	9649	9687	9711	9727
${\hat{\theta }}_y$	770	525	351	313	289	273

Table 5. Table 5: Useless auxiliary variable (Case 2): selection frequencies of θ ^ b subscript ^ 𝜃 𝑏 \hat{\theta}_{b} and θ ^ y subscript ^ 𝜃 𝑦 \hat{\theta}_{y} .

$n$	100	200	500	1000	2000	5000
${\hat{\theta }}_b$	1508	212	1	0	0	0
${\hat{\theta }}_y$	8492	9788	9999	10000	10000	10000

Table 6. Table 6: Useful auxiliary variable (Case 1): estimated risk functions of θ ^ b subscript ^ 𝜃 𝑏 \hat{\theta}_{b} , θ ^ y subscript ^ 𝜃 𝑦 \hat{\theta}_{y} and θ ^ b ​ e ​ s ​ t subscript ^ 𝜃 𝑏 𝑒 𝑠 𝑡 \hat{\theta}_{best} , and their standard errors in parenthesis

$n$	100	200	500	1000	2000	5000
$2n\left\{{ℛ}_x\left({\hat{\theta }}_b\right)-{ℒ}_x\left({\theta }_0\right)\right\}$	4.229	4.079	4.051	4.039	4.029	4.033
	(0.032)	(0.030)	(0.029)	(0.028)	(0.029)	(0.028)
$2n\left\{{ℛ}_x\left({\hat{\theta }}_y\right)-{ℒ}_x\left({\theta }_0\right)\right\}$	7.831	7.412	7.326	7.247	7.211	7.266
	(0.078)	(0.061)	(0.058)	(0.058)	(0.058)	(0.058)
$2n\left\{{ℛ}_x\left({\hat{\theta }}_{best}\right)-{ℒ}_x\left({\theta }_0\right)\right\}$	5.109	4.741	4.501	4.491	4.479	4.454
	(0.052)	(0.045)	(0.041)	(0.042)	(0.042)	(0.041)

Table 7. Table 7: Useless auxiliary variable (Case 2): estimated risk functions of θ ^ b subscript ^ 𝜃 𝑏 \hat{\theta}_{b} , θ ^ y subscript ^ 𝜃 𝑦 \hat{\theta}_{y} and θ ^ b ​ e ​ s ​ t subscript ^ 𝜃 𝑏 𝑒 𝑠 𝑡 \hat{\theta}_{best} , and their standard errors in parenthesis

$n$	100	200	500	1000	2000	5000
$2n\left\{{ℛ}_x\left({\hat{\theta }}_b\right)-{ℒ}_x\left({\theta }_0\right)\right\}$	105.527	214.659	543.685	1091.105	2182.647	5452.623
	(0.111)	(0.167)	(0.301)	(0.474)	(0.723)	(1.151)
$2n\left\{{ℛ}_x\left({\hat{\theta }}_y\right)-{ℒ}_x\left({\theta }_0\right)\right\}$	7.831	7.412	7.326	7.247	7.211	7.266
	(0.078)	(0.061)	(0.058)	(0.058)	(0.058)	(0.058)
$2n\left\{{ℛ}_x\left({\hat{\theta }}_{best}\right)-{ℒ}_x\left({\theta }_0\right)\right\}$	22.064	11.555	7.375	7.247	7.211	7.266
	(0.358)	(0.304)	(0.079)	(0.058)	(0.058)	(0.058)

Table 8. Table 8: Experiment average of n t ​ e ​ { ℒ ​ ( θ ^ y ) − ℒ x ​ ( θ ^ b ​ e ​ s ​ t ) } subscript 𝑛 𝑡 𝑒 ℒ subscript ^ 𝜃 𝑦 subscript ℒ 𝑥 subscript ^ 𝜃 𝑏 𝑒 𝑠 𝑡 n_{te}\{\mathcal{L}(\hat{\theta}_{y})-\mathcal{L}_{x}(\hat{\theta}_{best})\} for each case of Y = V ℓ 𝑌 subscript 𝑉 ℓ Y=V_{\ell} , ℓ = 1 , … , 13 ℓ 1 … 13 \ell=1,\ldots,13 . Standard errors are in parenthesis.

$Y$	$V_1$	$V_2$	$V_3$	$V_4$	$V_5$	$V_6$	$V_7$
$n_{te}\left\{{ℒ}_x\left({\hat{\theta }}_y\right)-{ℒ}_x\left({\hat{\theta }}_{best}\right)\right\}$	0.13	-0.14	89.71	46.24	-1.76	3.34	76.54
	(0.08)	(0.12)	(3.82)	(4.17)	(2.52)	(1.34)	(6.09)
$Y$	$V_8$	$V_9$	$V_{10}$	$V_{11}$	$V_{12}$	$V_{13}$
$n_{te}\left\{{ℒ}_x\left({\hat{\theta }}_y\right)-{ℒ}_x\left({\hat{\theta }}_{best}\right)\right\}$	13.91	39.45	1.72	111.24	15.48	0.23
	(2.21)	(3.12)	(0.29)	(8.46)	(2.11)	(0.09)

Equations

$${\hat{\theta }}_y$$

$$\Theta =\left\{\theta |\left(\begin{array}{c}\theta \\ \phi \end{array}\right)\in \mathcal{B}\right\}.$$

$${\hat{\beta }}_b$$

$$p_c(y,z,a;\beta )=p_x(y,z;\theta )p_{a|x}(a|y,z;\beta ).$$

$$D_x(q_x;p_x(\theta ))=-\int q_x(x)\log p_x(x;\theta )dx+\int q_x(x)\log q_x(x)dx.$$

$${\mathcal{L}}_x(\theta )=-\int q_x(x)\log p_x(x;\theta )dx.$$

$${\mathcal{R}}_x(\hat{\theta })=E[{\mathcal{L}}_x(\hat{\theta })],$$

$$q_{y|z}(y|z)$$

$$q_z(z)$$

$$\begin{array}{rl}{p}_{b|z}(y,a|z;\beta )& =z{N}_2((y,a{)}^{\top };{\mu }_1,\Sigma )+(1-z)N_2((y,a{)}^{\top };{\mu }_2,\Sigma ),\\ p_z(z;\theta )& ={\pi }_{1}z+(1-{\pi }_1)(1-z),\end{array}$$

$${\mu }_1=\left(\begin{array}{c}{\mu }_{1y}\\ {\mu }_{1a}\end{array}\right),{\mu }_2=\left(\begin{array}{c}{\mu }_{2y}\\ {\mu }_{2a}\end{array}\right),\Sigma =\left(\begin{array}{cc}{\sigma }_y^2& {\sigma }_{ya}\\ {\sigma }_{ya}& {\sigma }_a^2\end{array}\right).$$

$$\bar{\beta }=\left(\begin{array}{c}\bar{\theta }\\ \bar{\phi }\end{array}\right)=\arg \underset{\beta \in \mathcal{B}}{\max }\int q_b(b)\log p_b(b;\beta )db.$$

$$\sqrt{n}({\hat{\beta }}_b-\bar{\beta })\stackrel{d}{\to }{N}_{d+f}(0,I_b^{-1}{J}_b{I}_b^{-1}),$$

$$I_b=-E[{\nabla }^2\log p_b(b;\bar{\beta })],J_b=E[\nabla \log p_b(b;\bar{\beta }){\nabla }^{\top }\log p_b(b;\bar{\beta })].$$

$$I_x=-E[{\nabla }^2\log p_x(x;\bar{\theta })],I_y=-E[{\nabla }^2\log p_y(y;\bar{\theta })],I_{z|y}=-E[{\nabla }^2\log p_{z|y}(z|y;\bar{\theta })]=I_x-I_y,$$

$${\mathcal{R}}_x({\hat{\theta }}_b)$$

$$q_{z|y}(z|y)=p_{z|y}(z|y;\bar{\theta }).$$

$${\bar{\theta }}_y=\arg \underset{\theta \in \Theta }{\max }\int q_y(y)\log p_y(y;\theta )dy.$$

$$E[{\ell }_y({\hat{\theta }}_b)]$$

$${\mathcal{R}}_x({\hat{\theta }}_b)$$

$${\widehat{\mathrm{risk}}}_{x;b}$$

$$E[{\widehat{\mathrm{risk}}}_{x;b}]=2n{\mathcal{R}}_x({\hat{\theta }}_b)+2nC(q_x)+o(1).$$

$$\mathrm{TIC}=-2n{\ell }_y({\hat{\theta }}_y)+2\mathrm{tr}(I_y^{-1}{J}_y).$$

$$J_b=I_b,K_{b,y}=I_y,$$

$${\mathrm{AIC}}_{x;b}$$

$$E[{\mathrm{AIC}}_{x;b}]=2n{\mathcal{R}}_x({\hat{\theta }}_b)+2nC(q_x)+o(1).$$

$${\mathrm{AIC}}_{x;y}=-2n{\ell }_y({\hat{\theta }}_y)+\mathrm{tr}(I_x{I}_y^{-1})+d.$$

$${\mathrm{AIC}}_{y;b}$$

$${\mathrm{AIC}}_{y;y}=-2n{\ell }_y({\hat{\theta }}_y)+2d.$$

$${\mathrm{AIC}}_{x;b}-{\mathrm{AIC}}_{y;b}=\mathrm{tr}(I_x{I}_b^{-1})-\mathrm{tr}(I_y{I}_b^{-1})=\mathrm{tr}(I_{z|y}{I}_b^{-1})\ge 0,$$

$${\mathcal{L}}_x(\theta )$$

Full text

Akaike information criterion,

Auxiliary variables,

Fisher information matrix,

Incomplete data,

Kullback-Leibler divergence,

Misspecification,

Takeuchi information criterion,

and t1The research was supported in part by JSPS KAKENHI Grant (17K12650) and “Funds for the Development of Human Resources in Science and Technology” under MEXT, through the “Home for Innovative Researchers and Academic Knowledge Users (HIRAKU)” consortium. t2The research was supported in part by JSPS KAKENHI Grant (16H02789).

1 Introduction

Auxiliary variables are often observed along with primary variables. Here, the primary variables are random variables of interest, and our purpose is to estimate their predictive distribution, i.e., a probability distribution of the primary variables in future test data, while the auxiliary variables are random variables that are observed in training data but not included in the primary variables. We assume that the auxiliary variables are not observed in the test data, or we do not use them even if they are observed in the test data. When the auxiliary variables have a close relation with the primary variables, we expect to improve the accuracy of predictive distribution of the primary variables by considering a joint modeling of the primary and auxiliary variables.

The notion of auxiliary variables has been considered in statistics and machine learning literature. For example, the “curds and whey” method (Breiman and Friedman, 1997) and the “coaching variables” method (Tibshirani and Hinton, 1998) are based on a similar idea for improving prediction accuracy of primary variables by using auxiliary variables. In multitask learning, Caruana (1997) improved generalization accuracy of a main task by exploiting extra tasks. Auxiliary variables are also considered in incomplete data analysis, i.e., a part of primary variables are not observed; Mercatanti, Li and Mealli (2015) showed some theoretical results to make parameter estimation better by utilizing auxiliary variables in Gaussian mixture model (GMM).

Although auxiliary variables are expected to be useful for modeling primary variables, they can actually be harmful. As mentioned in Mercatanti, Li and Mealli (2015), using auxiliary variables may affect modeling result adversely because the number of parameters to be estimated increases and a candidate model of the auxiliary variables can be misspecified. Hence, it is important to select useful auxiliary variables. This is formulated as model selection by considering parametric models with auxiliary variables. In this paper, usefulness of auxiliary variables for estimating predictive distribution of primary variables is measured by a risk function based on the Kullback-Leibler (KL) divergence (Kullback and Leibler, 1951) that is often used for model selection. Because the KL risk function includes unknown parameters, we have to estimate it in actual use. Akaike Information Criterion (AIC) proposed by Akaike (1974) is one of the most famous criteria, which is known as an asymptotically unbiased estimator of the KL risk function. AIC is a good criterion from the perspective of prediction due to the asymptotic efficiency; see Shibata (1981, 1983). Takeuchi (1976) proposed a modified version of AIC, called Takeuchi Information Criterion (TIC), which relaxes an assumption for deriving AIC, that is, correct specification of candidate model. However, AIC and TIC are derived for primary variables without considering auxiliary variables in the setting of complete data analysis, and therefore, they are not suitable for auxiliary variable selection nor incomplete data analysis.

Incomplete data analysis is widely used in a broad range of statistical problems by regarding a part of primary variables as latent variables that are not observed. This setting also includes complete data analysis as a special case, where all the primary variables are observed. Information criteria for incomplete data analysis have been proposed in previous studies. Shimodaira (1994) developed an information criterion based on the KL divergence for complete data when the data are only partially observed. Cavanaugh and Shumway (1998) modified the first term of the information criterion of Shimodaira (1994) by the objective function of the EM algorithm (Dempster, Laird and Rubin, 1977). Recently, Shimodaira and Maeda (2018) proposed an information criterion, which is derived by mitigating a condition assumed in Shimodaira (1994) and Cavanaugh and Shumway (1998).

However, any of these previously proposed criteria are not derived by taking auxiliary variables into account. Thus, we propose a new information criterion by considering not only primary variable but also auxiliary variables in the setting of incomplete data analysis. The proposed criterion is a generalization of AIC, TIC and the criterion of Shimodaira and Maeda (2018). To the best of our knowledge, this is the first attempt to derive an information criterion by considering auxiliary variables. Moreover, we show an asymptotic equivalence between the proposed criterion and a variant of leave-one-out cross validation (LOOCV); this result is a generalization of the relationship between TIC and LOOCV (Stone, 1977).

Note that “auxiliary variables” may also be used in other contexts in literature. For example, Ibrahim, Lipsitz and Horton (2001) considered to use auxiliary variables in missing data analysis, which is similar to our usage in the sense that auxiliary variables are highly correlated with missing data. However, they use the auxiliary variables in order to avoid specifying a missing data mechanism; this goal is different from ours, because no missing data mechanism is considered in our study.

The reminder of this paper is organized as follows. Notations as well as a setting of this paper are introduced in Section 2. Illustrative examples of useful and useless auxiliary variables are given in Section 3. The information criterion for selecting useful auxiliary variables in incomplete data analysis is derived in Section 4, and the asymptotic equivalence between the proposed criterion and a variant of LOOCV is shown in Section 5. Performance of our method is examined via a simulation study and a real data analysis in Sections 6 and 7, respectively. Finally, we conclude this paper in Section 8. All proofs are shown in Appendix.

2 Preliminaries

2.1 Incomplete data analysis for primary variables

First we explain a setting of incomplete data analysis for primary variables in accordance with Shimodaira and Maeda (2018). Let $X$ denote a vector of primary variables, which consists of two parts as $X=(Y,Z)$, where $Y$ denotes the observed part and $Z$ denotes the unobserved latent part. This setting reduces to complete data analysis of $X=Y$ when $Z$ is empty. We write the true density function of $X$ as $q_{x}(x)=q_{x}(y,z)$ and a candidate parametric model of the true density as $p_{x}(x;\theta)=p_{x}(y,z;\theta)$, where $\theta\in\Theta\subset\mathbb{R}^{d}$ is an unknown parameter vector and $\Theta$ is its parameter space. We assume that $x=(y,z)\in\mathcal{Y}\times\mathcal{Z}$ for all density functions, where $\mathcal{Y}$ and $\mathcal{Z}$ are domains of $Y$ and $Z$, respectively. Thus the marginal densities of the observed part $Y$ are obtained by $q_{y}(y)=\int q_{x}(y,z)dz$ and $p_{y}(y;\theta)=\int p_{x}(y,z;\theta)dz$. For denoting densities, we will omit random variables such as $q_{y}$ and $p_{y}(\theta)$. We assume that $\theta$ is identifiable with respect to $p_{y}(\theta)$.

In this paper, we consider only a simple setting of i.i.d. random variables of sample size $n$. Let $x_{i}=(y_{i},z_{i})$, $i=1,\ldots,n$, be independent realizations of $X$, where we only observe $y_{1},\ldots,y_{n}$ and we cannot see the values of $z_{1},\ldots,z_{n}$. We estimate $\theta$ from the observed training data $y_{1},\ldots,y_{n}$. Then the maximum likelihood estimate (MLE) of $\theta$ is given by

[TABLE]

where $\ell_{y}(\theta)$ denotes the log-likelihood function (divided by $n$) of $\theta$ with respect to $y_{1},\ldots,y_{n}$.

If we were only interested in $Y$, we would consider the plug-in predictive distribution $p_{y}(\hat{\theta}_{y})$ by substituting $\hat{\theta}_{y}$ into $p_{y}(\theta)$. However, we are interested in the whole primary variable $X=(Y,Z)$ and its density $q_{x}$. We thus consider $p_{x}(\hat{\theta}_{y})$ by substituting $\hat{\theta}_{y}$ into $p_{x}(\theta)$, and evaluate the MLE by comparing $p_{x}(\hat{\theta}_{y})$ with $q_{x}$. For this purpose, Shimodaira and Maeda (2018) derived an information criterion as an asymptotically unbiased estimator of the KL risk function which measures how well $p_{x}(\hat{\theta}_{y})$ approximates $q_{x}$.

2.2 Statistical analysis with auxiliary variables

Next, we extend the setting to incomplete data analysis with auxiliary variables. Let $A$ denote a vector of auxiliary variables. In addition to $Y$, we observe $A$ in the training data, but we are not interested in $A$. For convenience, we introduce a vector of observable variables $B=(Y,A)$ and a vector of all variables $C=(Y,Z,A)$ as summarized in Table 1. Now $c_{i}=(y_{i},z_{i},a_{i})$, $i=1,\ldots,n$, are independent realizations of $C$, and we estimate $\theta$ from the observed training data $b_{i}=(y_{i},a_{i})$, $i=1,\ldots,n$. Let $\hat{\theta}_{b}$ be the MLE of $\theta$ by using $A$ in addition to $Y$. Since we are only interested in the primary variables, we consider the plug-in predictive distribution $p_{x}(\hat{\theta}_{b})$ by substituting $\hat{\theta}_{b}$ into $p_{x}(\theta)$, and evaluate the MLE by comparing $p_{x}(\hat{\theta}_{b})$ with $q_{x}$.

In order to define the MLE $\hat{\theta}_{b}$, let us clarify a candidate parametric model with auxiliary variables. We write the true density function of $C$ as $q_{c}(c)=q_{c}(y,z,a)$ and a candidate parametric model of the true density as $p_{c}(c;\beta)=p_{c}(y,z,a;\beta)$, where $\beta=(\theta^{\top},\varphi^{\top})^{\top}\in\mathcal{B}\subset\mathbb{R}^{d+f}$ is an unknown parameter vector with nuisance parameter $\varphi\in\mathbb{R}^{f}$ and $\mathcal{B}$ is its parameter space. We assume that $c=(y,z,a)\in\mathcal{Y}\times\mathcal{Z}\times\mathcal{A}$ for all density functions, where $\mathcal{A}$ is the domain of $A$. We also assume that $\beta$ is identifiable with respect to $p_{b}(y,a;\beta)=\int p_{c}(y,z,a;\beta)dz$. Let us redefine $p_{x}(\theta)$ as $p_{x}(y,z;\theta)=\int p_{c}(y,z,a;\beta)da$ and the parameter space of $\theta$ as

[TABLE]

Then, $\hat{\theta}_{b}$ is obtained from the MLE of $\beta$ given by

[TABLE]

where $\ell_{b}(\beta)$ denotes the log-likelihood function (divided by $n$) of $\beta$ with respect to $b_{1},\ldots,b_{n}$.

Finally, we introduce a general notation for density functions. For a random variable, say $R$, we write the true density function as $q_{r}(r)$ and a candidate parametric model of $q_{r}$ as $p_{r}(r;\theta)$ or $p_{r}(r;\beta)$. For random variables $R$ and $S$, we write the true conditional density function of $R$ given $S=s$ as $q_{r|s}(r|s)$ and its corresponding model as $p_{r|s}(r|s;\theta)$ or $p_{r|s}(r|s;\beta)$. For example, a candidate model of $C$ can be decomposed as

[TABLE]

2.3 Comparing the two estimators

We have thus far obtained the two MLEs of $\theta$, namely $\hat{\theta}_{y}$ and $\hat{\theta}_{b}$, and their corresponding predictive distributions $p_{x}(\hat{\theta}_{y})$ and $p_{x}(\hat{\theta}_{b})$, respectively. We would like to determine which of the two predictive distributions approximates $q_{x}$ better than the other. The approximation error of $p_{x}(\theta)$ is measured by the KL divergence from $q_{x}$ to $p_{x}(\theta)$ defined as

[TABLE]

Since the last term on the right hand side does not depend on $p_{x}(\theta)$, we ignore it for computing the loss function of $p_{x}(\theta)$ defined by

[TABLE]

Let $\hat{\theta}$ be an estimator of $\theta$. The risk (or expected loss) function of $p_{x}(\hat{\theta})$ is defined by

[TABLE]

where we take the expectation by considering $\hat{\theta}$ as a random variable. Note that $\hat{\theta}$ in the notation of $\mathcal{R}_{x}(\hat{\theta})$ indicates the procedure for computing $\hat{\theta}$ instead of a particular value of $\hat{\theta}$. $\mathcal{R}_{x}(\hat{\theta})$ measures how well $p_{x}(\hat{\theta})$ approximates $q_{x}$ on average in the long run.

For comparing the two MLEs, we define $\mathcal{R}_{x}(\hat{\theta}_{y})$ and $\mathcal{R}_{x}(\hat{\theta}_{b})$ by considering that $\hat{\theta}_{y}$ and $\hat{\theta}_{b}$ are functions of independent random variables $Y_{1},\ldots,Y_{n}$ and $B_{1},\ldots,B_{n}$, respectively, where $B_{i}=(Y_{i},A_{i})$ has the same distribution as $B$ for all $i=1,\ldots,n$. $\hat{\theta}_{b}$ is better than $\hat{\theta}_{y}$ when $\mathcal{R}_{x}(\hat{\theta}_{b})<\mathcal{R}_{x}(\hat{\theta}_{y})$, that is, the auxiliary variable $A$ helps the statistical inference on $q_{x}$. On the other hand, $A$ is harmful when $\mathcal{R}_{x}(\hat{\theta}_{b})>\mathcal{R}_{x}(\hat{\theta}_{y})$. Although we focus only on comparison between $Y$ and $B=(Y,A)$ in this paper, if there are more than two auxiliary variables (and their combinations) $A_{1},A_{2},\ldots$, then we may compare $\mathcal{R}_{x}(\hat{\theta}_{(y,a_{1})}),\mathcal{R}_{x}(\hat{\theta}_{(y,a_{2})}),\ldots$, to determine good auxiliary variables. Of course, the risk functions cannot be calculated in reality because they depend on the unknown true distribution. Thus, we derive a new information criterion as an estimator of the risk function in our setting. Since an asymptotically unbiased estimator of $\mathcal{R}_{x}(\hat{\theta}_{y})$ has been already derived in Shimodaira and Maeda (2018), we will only derive an asymptotically unbiased estimator of $\mathcal{R}_{x}(\hat{\theta}_{b})$.

3 An illustrative example with auxiliary variables

3.1 Model setting

In this section, we demonstrate parameter estimation by using auxiliary variables in Gaussian mixture model (GMM), which can be formulated in incomplete data analysis. Let us consider two-component GMM; observed values are generated from one of two Gaussian distributions, where the assigned labels are missing. The observed data and missing labels are realizations of $Y$ and $Z$, respectively. We estimate a predictive distribution of $X=(Y,Z)$ from the observation of $Y$, and we attempt improving it by utilizing $A$ in addition to $Y$. The true density function of primary variables $X=(Y,Z)\in\mathbb{R}\times\{0,1\}$ is given as

[TABLE]

where $N(\cdot;\mu,\sigma^{2})$ denotes the density function of $N(\mu,\sigma^{2})$, i.e., the normal distribution with mean $\mu$ and variance $\sigma^{2}$. We consider the following two cases for the true conditional distribution of auxiliary variable $A$ given $X=x$:

• • Case 1:

$q_{a|x}(a|y,z)=q_{a|z}(a|z)=zN(a;1.8,0.49)+(1-z)N(a;-1.8,0.49)$.

• Case 2:

$q_{a|x}(a|y,z)=q_{a}(a)=0.6N(a;1.8,0.49)+0.4N(a;-1.8,0.49)$.

The random variables $X$ and $A$ are not independent in Case 1 whereas they are independent in Case 2. Hence, $A$ will contribute to estimating $\theta$ in Case 1. On the other hand, in Case 2, $A$ must not be useful, and $A$ becomes just noise if we estimate $\theta$ from $Y$ and $A$.

In both cases, we use the following two-component GMM as a candidate model of $q_{c}$:

[TABLE]

where $N_{2}(\cdot;\mu_{i},\Sigma)$ denotes the density function of bivariate normal distribution $N_{2}(\mu_{i},\Sigma)$, $i=1,2$, and the parameters are

[TABLE]

Therefore, $\beta=(\theta^{\top},\varphi^{\top})^{\top}$, $\theta=(\pi_{1},\mu_{1y},\mu_{2y},\sigma_{y}^{2})^{\top}$ and $\varphi=(\mu_{1a},\mu_{2a},\sigma_{a}^{2},\sigma_{ya})^{\top}$. The true parameters of $\theta$ and $\varphi$ for Case 1 are given by $\theta_{0}=(0.6,-1.2,1.2,0.7)^{\top}$ and $\varphi_{0}=(1.8,-1.8,0.49,0)^{\top}$, respectively. By considering the joint density function $p_{c}(y,z,a;\beta)=p_{b|z}(y,a|z;\beta)p_{z}(z;\theta)$, this candidate model correctly specifies the true density function $q_{c}(y,z,a)=q_{a|x}(a|y,z)q_{y|z}(y|z)q_{z}(z)$ in Case 1. On the other hand, the model is misspecified for Case 2, and we cannot think of the true parameters.

3.2 Estimation results

For illustrating the impact of auxiliary variables on parameter estimation in each case, we generated a typical dataset $c_{1},\ldots,c_{n}$ with sample size $n=100$ from $q_{c}$, which is actually picked from 10,000 datasets generated in the simulation study of Section 6, and details of how to select the typical dataset are also shown there. For each case, we computed the three MLEs $\hat{\theta}_{y}$, $\hat{\theta}_{b}$ and $\hat{\theta}_{x}$, where $\hat{\theta}_{x}$ is the MLE of $\theta$ calculated by using complete data $x_{1},\ldots,x_{n}$ as if labels $z_{1},\ldots,z_{n}$ were available.

The result of Case 1 is shown in Figure 1, where $A$ is beneficial for estimating $\theta$. In the left panel, the two clusters are well separated, which makes parameter estimation stable. The estimated $p_{b}(\hat{\beta}_{b})$ captures the structure of the two clusters corresponding to the label $z_{i}=0$ and $z_{i}=1$, showing that $p_{c}(\hat{\beta}_{b})$ is estimated reasonably well, and thus $p_{x}(\hat{\theta}_{b})$ is a good approximation of $q_{x}$. Looking at the right panel, we also observe that $p_{y}(\hat{\theta}_{b})$ is better than $p_{y}(\hat{\theta}_{y})$ for approximating $p_{y}(\hat{\theta}_{x})$, suggesting that the auxiliary variable is useful for recovering the lost information of missing data. In fact, the three MLEs are calculated as follows: $\hat{\theta}_{y}=(0.671,-1.143,1.324,0.678)^{\top}$, $\hat{\theta}_{b}=(0.613,-1.228,1.093,0.744)^{\top}$ and $\hat{\theta}_{x}=(0.620,-1.233,1.141,0.695)^{\top}$. By comparing $\|\hat{\theta}_{b}-\hat{\theta}_{x}\|=0.069$ with $\|\hat{\theta}_{y}-\hat{\theta}_{x}\|=0.212$, we can see that $\hat{\theta}_{b}$ is better than $\hat{\theta}_{y}$ for predicting $\hat{\theta}_{x}$ without looking at the latent variable. All these observations indicate that the parameter estimation of $\theta$ is improved by using $A$ in Case 1.

The result of Case 2 is shown in Figure 2, where $A$ is harmful for estimating $\theta$. For fair comparison, exactly the same values of $\{(y_{i},z_{i})\}_{i=1}^{100}$ are used in both cases. Thus, $\hat{\theta}_{y}$ and $\hat{\theta}_{x}$ have the same values as in Case 1 whereas $\hat{\theta}_{b}$ has a different value as $\hat{\theta}_{b}=(0.581,-0.403,-0.232,2.015)^{\top}$. By comparing $\|\hat{\theta}_{b}-\hat{\theta}_{x}\|=2.078$ with $\|\hat{\theta}_{y}-\hat{\theta}_{x}\|=0.212$, we can see that $\hat{\theta}_{b}$ is worse than $\hat{\theta}_{y}$ for predicting $\hat{\theta}_{x}$. This is also seen in Figure 2. In the left panel, the estimated $p_{b}(\hat{\beta}_{b})$ captures some structure of the two clusters, but they do not correspond to the label $z_{i}=0$ and $z_{i}=1$. As a result, $p_{y}(\hat{\theta}_{b})$ becomes a very poor approximation of $p_{y}(\hat{\theta}_{x})$ in the right panel, indicating that the parameter estimation of $\theta$ is actually hindered by using $A$ in Case 2.

These examples suggest that usefulness of auxiliary variables depends strongly on the true distribution and a candidate model. Hence, it is important to select useful auxiliary variables from observed data.

4 Information criterion

4.1 Asymptotic expansion of the risk function

In this section, we derive a new information criterion as an asymptotically unbiased estimator of the risk function $\mathcal{R}_{x}(\hat{\theta}_{b})$ defined in (3). We start from a general framework of misspecification, i.e., without assuming that candidate models are correctly specified, and later we give specific assumptions. Let $\bar{\beta}$ be the optimal parameter value with respect to the KL divergence from $q_{b}$ to $p_{b}(\beta)$, that is,

[TABLE]

If the candidate model is correctly specified, i.e., there exists $\beta_{0}=(\theta_{0}^{\top},\varphi_{0}^{\top})^{\top}$ such that $q_{b}=p_{b}(\beta_{0})$, then $\bar{\beta}=\beta_{0}$ as well as $\bar{\theta}=\theta_{0}$.

In this paper, we assume the regularity conditions A1 to A6 of White (1982) for $q_{b}$ and $p_{b}(\beta)$ so that the MLE $\hat{\beta}_{b}$ has consistency and asymptotic normality. In particular, $\bar{\beta}$ is determined uniquely (i.e., identifiable) and is interior to $\mathcal{B}$. We assume that $I_{b}$ and $J_{b}$ defined below are nonsingular in the neighbourhood of $\bar{\beta}$. Then White (1982) showed the asymptotic normality as $n\to\infty$,

[TABLE]

where $I_{b}$ and $J_{b}$ are $(d+f)\times(d+f)$ matrices defined by using $\nabla=\partial/\partial\beta$, $\nabla^{\top}=\partial/\partial\beta^{\top}$ and $\nabla^{2}=\partial^{2}/\partial\beta\partial\beta^{\top}$ as

[TABLE]

Note that we write derivatives by abbreviated forms, e.g., $\nabla^{2}\log{p_{b}}(b;\bar{\beta})$ means $\nabla^{2}\log{p_{b}}(b;\beta)|_{\beta=\bar{\beta}}$ and so on. In addition, we allow interchange of integrals and derivatives rather formally when working with models, although we actually need conditions for the models such as White (1982). Moreover, the condition A7 of White (1982) is assumed in order to establish $I_{b}=J_{b}$ when considering a situation that the candidate model is correctly specified. We assume the above conditions throughout the paper without explicitly stated.

Let us define three $(d+f)\times(d+f)$ matrices as

[TABLE]

which will be used in the lemmas below. Since the derivatives of $\log{p_{x}}(x;\theta)$ and $\log{p_{y}}(y;\theta)$ with respect to $\varphi$ is zero, the matrices become singular when $f>0$, but this is not a problem in our calculation. The following lemma shows that the dominant term of $\mathcal{R}_{x}(\hat{\theta}_{b})$ is $\mathcal{L}_{x}(\bar{\theta})$ and the remainder terms are of order $O(n^{-1})$, by noting that $\nabla^{\top}\mathcal{L}_{x}(\bar{\theta})=O(1)$ and $E[\hat{\beta}_{b}-\bar{\beta}]=O(n^{-1})$ in general. The proof is given in Appendix A.1.

Lemma 1**.**

The risk function $\mathcal{R}_{x}(\hat{\theta}_{b})$ is expanded asymptotically as

[TABLE]

Just as a remark, the term $\nabla^{\top}\mathcal{L}_{x}(\bar{\theta})E[\hat{\beta}_{b}-\bar{\beta}]=O(n^{-1})$ above does not appear in the derivation of AIC or TIC, where $B=X$ and thus $\nabla^{\top}\mathcal{L}_{x}(\bar{\theta})=0$. This term appears when the loss function for evaluation and that for estimation differ, for example, in the derivation of the information criterion under covariate shift; see $K_{w}^{[1]\top}b_{w}$ in eq. (4.1) of Shimodaira (2000).

4.2 Estimating the risk function

For deriving an estimator of $\mathcal{R}_{x}(\hat{\theta}_{b})$, we introduce an additional condition. Let us assume that the candidate model is correctly specified for the latent part as

[TABLE]

This is the same condition as eq. (14) of Shimodaira and Maeda (2018) except that $\bar{\theta}$ is replaced by

[TABLE]

Since $Z$ is missing completely in our setting, we need such a condition to proceed further. Although any method cannot detect misspecification of $p_{z|y}$ if $p_{b}$ is correctly specified, it is often the case that misspecification of $p_{z|y}$ leads to that of $p_{b}$, and thus it is detected indirectly as in Case 2 of Section 3.

Note that the symbol of $\bar{\theta}$ in our notation should have been $\bar{\theta}_{b}$, although we used $\bar{\theta}$ for simplicity, and there is also $\bar{\theta}_{x}$ defined similarly from $p_{x}(x;\theta)$. They all differ each other with differences of order $O(1)$ in general, but $\bar{\theta}=\bar{\theta}_{y}=\bar{\theta}_{x}=\theta_{0}$ when $p_{c}(\beta)$ is correctly specified as $q_{c}=p_{c}(\beta_{0})$.

Now we give the asymptotic expansion of $E[\ell_{y}(\hat{\theta}_{b})]$, which shows that $-\ell_{y}(\hat{\theta}_{b})$ can be used as an estimator of $\mathcal{L}_{x}(\bar{\theta})$ but the asymptotic bias is of order $O(n^{-1})$.

Lemma 2**.**

Assume the condition (6). Then, the expectation of the estimated log-likelihood $\ell_{y}(\hat{\theta}_{b})$ can be expanded as

[TABLE]

where $K_{b,y}=E[\nabla\log{p_{b}}(\bar{\beta})\nabla^{\top}\log{p_{y}}(\bar{\theta})]$ and $C(q_{x})=\int q_{x}(x)\log{q_{z|y}}(z|y)dx$.

The proof of Lemma 2 is given in Appendix A.2. By eliminating $\mathcal{L}_{x}(\bar{\theta})$ from the two expressions in Lemma 1 and Lemma 2, and rearranging the formula, we get the following lemma, which plays a central role in deriving our information criterion.

Lemma 3**.**

Assume the condition (6). Then, an expansion of the risk function $\mathcal{R}_{x}(\hat{\theta}_{b})$ is given by

[TABLE]

We can ignore $C(q_{x})$ for model selection, because it is a constant term which does not depend on the candidate model. Thus, finally, we define an information criterion from the right hand side of (7). The following theorem is an immediate consequence of Lemma 3.

Theorem 4**.**

Assume the condition (6). Let us define an information criterion as

[TABLE]

Then this criterion is an asymptotically unbiased estimator of $2n\mathcal{R}_{x}(\hat{\theta}_{b})$ by ignoring the constant term $C(q_{x})$.

[TABLE]

Note that the subscript of $\widehat{{\rm risk}}_{x;b}$, $x;b$ is defined in accordance with Shimodaira and Maeda (2018); thus the former $x$ and the latter $b$ mean random variables used in evaluation and estimation, respectively. This criterion is an extension of TIC because when $X=B=Y$, $\widehat{{\rm risk}}_{x;b}$ coincides with TIC of Takeuchi (1976) defined as follows:

[TABLE]

4.3 Akaike information criteria for auxiliary variable selection

In actual use, $\widehat{{\rm risk}}_{x;b}$ may have a too complicated form. Thus, we derive a simpler information criterion by assuming the correctness of the candidate model like as AIC.

Theorem 5**.**

Suppose $p_{c}(\beta)$ is correctly specified so that $q_{c}=p_{c}(\beta_{0})$ for some $\beta_{0}\in\mathcal{B}$. Then, we have

[TABLE]

and thus $\widehat{{\rm risk}}_{x;b}$ is rewritten as

[TABLE]

This criterion is an asymptotically unbiased estimator of $2n\mathcal{R}_{x}(\hat{\theta}_{b})$ by ignoring the constant term $C(q_{x})$.

[TABLE]

The proof is given in Appendix A.3. $I_{x}$, $I_{y}$ and $I_{b}$ are replaced by their consistent estimators in practical situations.

The newly obtained criterion ${\rm AIC}_{x;b}$ is a generalization of AIC and some of its variants. If $\theta$ is estimated by $\hat{\theta}_{y}$ instead of $\hat{\theta}_{b}$, we simply let $B=Y$ in the expression of ${\rm AIC}_{x;b}$ so that we get ${\rm AIC}_{x;y}$ proposed by Shimodaira and Maeda (2018):

[TABLE]

Note that if $B=Y$, $I_{y}$ is not singular because $\beta=\theta$. On the other hand, if there is no latent part, we simply let $X=Y$ in the expression of ${\rm AIC}_{x;b}$ so that we get

[TABLE]

This can be used to select useful auxiliary variables in complete data analysis. Moreover, if $X=Y=B$, ${\rm AIC}_{x;b}$ reduces to the original AIC proposed by Akaike (1974):

[TABLE]

It is worth mentioning that ${\rm tr}(I_{z|y}I_{b}^{-1})$ is interpreted as the additional penalty for the latent part:

[TABLE]

which is also mentioned in eq. (1) of Shimodaira and Maeda (2018) for the case of $B=Y$.

4.4 The illustrative example (cont.)

Let us return to the problem of determining whether to use the auxiliary variables or not, that is, comparison between $p_{x}(\hat{\theta}_{b})$ and $p_{x}(\hat{\theta}_{y})$. By comparing ${\rm AIC}_{x;b}$ with ${\rm AIC}_{x;y}$, we can determine whether the vector of auxiliary variables $A$ is useful or useless. Thus, only when ${\rm AIC}_{x;b}<{\rm AIC}_{x;y}$, we conclude that $A$ is useful in order to estimate $\theta$ for predicting $X$.

Let us apply this procedure to the illustrative example in Section 3. The generalized AICs are computed for the two cases of the typical dataset, and the results are shown in Table 2. Looking at the value of ${\rm AIC}_{x;b}-{\rm AIC}_{x;y}$, it is negative for Case 1, concluding that the auxiliary variable is useful, and it is positive for Case 2, concluding that the auxiliary variable is useless. According to the AIC values, therefore, we use the auxiliary variable of Case 1, but do not use the auxiliary variable of Case 2. This decision agrees with the observations of Figures 1 and 2 in Section 3.2. Actually, the decision is correct, because the value of $\mathcal{R}_{x}(\hat{\theta}_{b})-\mathcal{R}_{x}(\hat{\theta}_{y})$ is negative for Case 1 and positive for Case 2 as will be seen in the simulation study of Section 6.2.

We can also argue the usefulness of the auxiliary variable for predicting $Y$ instead of $X$, that is, comparison between $p_{y}(\hat{\theta}_{b})$ and $p_{y}(\hat{\theta}_{y})$. By comparing ${\rm AIC}_{y;b}$ with ${\rm AIC}_{y;y}$, we can determine whether $A$ is useful or useless for predicting $Y$. Looking at the value of ${\rm AIC}_{y;b}-{\rm AIC}_{y;y}$ in Table 2, we make the same decision as that for $X$.

5 Leave-one-out cross validation

Variable selection by cross-validatory (CV) choice (Stone, 1974) is often applied to real data analysis due to its simplicity, although its computational burden is larger than that of information criteria; see Arlot and Celisse (2010) for a recent review of cross-validation methods. As shown in Stone (1977), leave-one-out cross validation (LOOCV) is asymptotically equivalent to TIC. Because LOOCV does not require calculation of the information matrices of TIC, LOOCV is easier to use than TIC. There are also some literature for improving LOOCV such as Yanagihara, Tonda and Matsumoto (2006), which gives a modification of LOOCV to reduce its bias by considering maximum weighted log-likelihood estimation. However, we focus on the result of Stone (1977) and extend it to our setting.

In incomplete data analysis, LOOCV cannot be directly used because the loss function with respect to the complete data includes latent variables. Thus, we transform the loss function as follows:

[TABLE]

where $g(y;\theta)=\log p_{y}(y;\theta)+f(y;\theta)$ and

[TABLE]

Note that $f(y;\theta)=0$ when $X=Y$. Using the function $g(y;\theta)$, we then obtain the following LOOCV estimator of the risk function $\mathcal{R}_{x}(\hat{\theta}_{b})$.

[TABLE]

where $\hat{\theta}_{b}^{(-i)}$ is the leave-out-out estimate of $\theta$ defined as

[TABLE]

We will show below in this section that $\mathcal{L}_{x}^{\rm cv}(\hat{\theta}_{b})$ is asymptotically equivalent to $\widehat{{\rm risk}}_{x;b}$. For implementing the LOOCV procedure with latent variables, however, we have to estimate $q_{z|y}(z|y)$ by $p_{z|y}(z|y,\hat{\theta}_{b})$ in $f(y;\theta)$. This introduces a bias to $\mathcal{L}_{x}^{\rm cv}(\hat{\theta}_{b})$, and hence, information criteria are preferable to the LOOCV in incomplete data analysis.

Let us show the asymptotic equivalence of $\mathcal{L}_{x}^{\rm cv}(\hat{\theta}_{b})$ and $\widehat{{\rm risk}}_{x;b}$ by assuming that we know the functional form of $f(y;\theta)$. Noting that $\hat{\beta}_{b}^{(-i)}$ is a critical point of $\ell_{b}(\beta)-\log p_{b}(b_{i};\beta)/n$, we have

[TABLE]

By applying Taylor expansion to $\nabla\ell_{b}(\beta)$ around $\beta=\hat{\beta}_{b}$, it follows from $\nabla\ell_{b}(\hat{\beta}_{b})=0$ that

[TABLE]

where $\tilde{\beta}_{b}^{i}$ lies between $\hat{\beta}_{b}^{(-i)}$ and $\hat{\beta}_{b}$. We can see from (14) that $\hat{\beta}_{b}^{(-i)}-\hat{\beta}_{b}=O_{p}(n^{-1})$. Next, we regard $g(y_{i};\theta)$ as a function of $\beta$ and apply Taylor expansion to it around $\beta=\hat{\beta}_{b}$. Therefore, $g(y_{i};\hat{\theta}_{b}^{(-i)})$ can be expressed as follows:

[TABLE]

where $\tilde{\theta}_{b}^{i}$ lies between $\hat{\theta}_{b}^{(-i)}$ and $\hat{\theta}_{b}$ ($\tilde{\theta}_{b}^{i}$ does not corresponding to $\tilde{\beta}_{b}^{i}$). Then we assume that

[TABLE]

By noting $\hat{\beta}_{b}^{(-i)}=\hat{\beta}_{b}+O_{p}(n^{-1})$, we have $\tilde{\beta}_{b}^{i}=\bar{\beta}+O_{p}(n^{-1/2})$ and $\tilde{\theta}_{b}^{i}=\bar{\theta}+O_{p}(n^{-1/2})$, and thus (16) holds at least formally. With the above setup, we show the following theorem. The proof is given in Appendix A.4.

Theorem 6**.**

Suppose the same assumptions of Theorem 4 and (16), we have

[TABLE]

Because the second term on the right-hand side of (17) does not depend on candidate models under the condition (6), this theorem implies that $\mathcal{L}_{x}^{\rm cv}(\hat{\theta}_{b})$ is asymptotically equivalent to $\widehat{{\rm risk}}_{x;b}$ except for the scaling and the constant term. However, someone may wonder why $f(y;\theta)$ is included in $g(y;\theta)$ for comparing models of $p(b;\beta)$. By assuming that $p_{z|y}(\theta)$ is correctly specified for $q_{z|y}$, $f(y;\bar{\theta})=\int q_{z|y}(z|y)\log q_{z|y}(z|y)dz$ does not depend on the model anymore, so we may simply exclude $f(y;\theta)$ from $g(y;\theta)$, leading to the loss $\mathcal{L}_{y}(\theta)$ instead. The reason of including $f(y;\theta)$ in $g(y;\theta)$ is explained as follows. $\mathcal{L}_{x}^{\rm cv}(\hat{\theta}_{b})$ as well as $\widehat{{\rm risk}}_{x;b}$ (and ${\rm AIC}_{x;b}$) includes the additional penalty for estimating $\hat{\theta}_{b}$ in $f(y;\hat{\theta}_{b})$, which depends on the candidate models even if $p_{z|y}(\theta)$ is correctly specified.

6 Experiments with simulated datasets

This section shows the usefulness of auxiliary variables and the proposed information criteria via a simulation study. The models illustrated in Section 3 are used for confirming the asymptotic unbiasedness of the information criterion and the validity of auxiliary variable selection.

6.1 Unbiasedness

At first, we confirm the asymptotic unbiasedness of ${\rm AIC}_{x;b}$ for estimating $2n\mathcal{R}_{x}(\hat{\theta}_{b})$ except for the constant term, $C(q_{x})$. The simulation setting is the same as Case 1 in Section 3, thus the data generating model is given by

[TABLE]

where $\mu_{10}=-\mu_{20}=(-1.2,1.8)^{\top}$ and $\Sigma_{0}={\rm diag}(0.7,0.49)$. We generated $T=10^{4}$ independent replicates of the dataset $\{(y_{i},z_{i},a_{i})\}_{i=1}^{n}$ from this model; in fact, we used $\{(y_{i},z_{i},a_{i,1})\}_{i=1}^{n}$ generated in Section 6.2. The candidate model is given by (4), which is correctly specified for the above data generating model. Because ${\rm AIC}_{x;b}$ is derived by ignoring $C(q_{x})$, we compare $E[{\rm AIC}_{x;b}-{\rm AIC}_{x;y}]$ with $2n\{\mathcal{R}_{x}(\hat{\theta}_{b})-\mathcal{R}_{x}(\hat{\theta}_{y})\}$. The computation of the expectation is approximated by the simulation average as

[TABLE]

where ${\rm AIC}_{x;b}^{(t)}$, ${\rm AIC}_{x;y}^{(t)}$, $\hat{\theta}_{b}^{(t)}$ and $\hat{\theta}_{y}^{(t)}$ are those computed for the $t$-th dataset ($t=1,\ldots,T$).

Here, we remark about calculation of the loss function $\mathcal{L}_{x}(\hat{\theta})$ in two-component GMM. Let $\hat{\theta}=(\hat{\pi}_{1},\hat{\mu}_{1},\hat{\mu}_{2},\hat{\sigma}^{2})^{\top}$ be an estimator of $\theta$. We expect that the components of GMM corresponding to $Z=1$ and $Z=0$ consists of $(\hat{\pi}_{1},\hat{\mu}_{1},\hat{\sigma}^{2})$ and $(1-\hat{\pi}_{1},\hat{\mu}_{2},\hat{\sigma}^{2})$, respectively. However, we cannot determine the assignment of the estimated parameters in reality, i.e., $(\hat{\pi}_{1},\hat{\mu}_{1},\hat{\sigma}^{2})$ and $(1-\hat{\pi}_{1},\hat{\mu}_{2},\hat{\sigma}^{2})$ may correspond to $Z=0$ and $Z=1$, respectively, because the labels $z_{1},\ldots,z_{n}$ are missing. The assignment is required to calculate $\mathcal{L}_{x}(\hat{\theta})$ whereas it is not used for $\mathcal{L}_{y}(\hat{\theta})$ and the proposed information criteria. Hence, in this paper, we define $\mathcal{L}_{x}(\hat{\theta})$ as the minimum value between $\mathcal{L}(\hat{\theta})$ and $\mathcal{L}(\hat{\theta}^{\prime})$, where $\hat{\theta}^{\prime}=(1-\hat{\pi}_{1},\hat{\mu}_{2},\hat{\mu}_{1},\hat{\sigma}^{2})^{\top}$.

Table 3 shows the result of the simulation for $n=100,200,500,1000,2000$ and $5000$. For all $n$, we observe that $E[{\rm AIC}_{x;b}-{\rm AIC}_{x;y}]$ is very close to $2n\{\mathcal{R}_{x}(\hat{\theta}_{b})-\mathcal{R}_{x}(\hat{\theta}_{y})\}$, indicating the unbiasedness of ${\rm AIC}_{x;b}$.

6.2 Auxiliary variable selection

Next, we demonstrate that the proposed AIC selects a useful auxiliary variable (Case 1), while it does not select a useless auxiliary variable (Case 2). In each case, we generated $T=10^{4}$ independent replicates of the dataset $\{(y_{i},z_{i},a_{i})\}_{i=1}^{n}$ from the model. Actually, the values of $\{(y_{i},z_{i})\}_{i=1}^{n}$ are shared in both cases, so we generated replicates of $\{(y_{i},z_{i},a_{i,1},a_{i,2})\}_{i=1}^{n}$, where $a_{i,1}$ and $a_{i,2}$ are auxiliary variables for Case 1 and Case 2, respectively. In each case, we compute ${\rm AIC}_{x;b}$ and ${\rm AIC}_{x;y}$, then we select $\hat{\theta}_{b}$ (i.e., selecting the auxiliary variable $A$) if ${\rm AIC}_{x;b}<{\rm AIC}_{x;y}$ and select $\hat{\theta}_{y}$ (i.e., not selecting the auxiliary variable $A$) otherwise. The selected estimator is denoted as $\hat{\theta}_{best}$. This experiment was repeated for $T=10^{4}$ times. Note that the typical dataset in Section 3 was picked from the generated datasets so that it has around the median value in each of $\mathcal{L}_{x}(\hat{\theta}_{b})-\mathcal{L}_{x}(\hat{\theta}_{y})$, $\mathcal{L}_{y}(\hat{\theta}_{b})-\mathcal{L}_{y}(\hat{\theta}_{y})$, ${\rm AIC}_{x;b}-{\rm AIC}_{x;y}$ and ${\rm AIC}_{y;b}-{\rm AIC}_{y;y}$ in both cases.

The selection frequencies are shown in Tables 4 and 5. We observe that, as expected, the useful auxiliary variable tends to be selected in Case 1, while the useless auxiliary variable tends to be not selected in Case 2.

For verifying the usefulness of the auxiliary variable in both cases, we computed the risk value $\mathcal{R}_{x}(\hat{\theta})$ for $\hat{\theta}=\hat{\theta}_{y}$, $\hat{\theta}_{b}$ and $\hat{\theta}_{best}$. They are approximated by the simulation average as

[TABLE]

The results are shown in Tables 6 and 7. For easier comparisons, the values are the differences from $\mathcal{L}_{x}(\theta_{0})$ with the true value $\theta_{0}$. For all $n$, we observe that, as expected, $\mathcal{R}_{x}(\hat{\theta}_{b})<\mathcal{R}_{x}(\hat{\theta}_{y})$ in Case 1, and $\mathcal{R}_{x}(\hat{\theta}_{b})>\mathcal{R}_{x}(\hat{\theta}_{y})$ in Case 2. In both cases, $\mathcal{R}_{x}(\hat{\theta}_{best})$ is close to $\min\{\mathcal{R}_{x}(\hat{\theta}_{b}),\mathcal{R}_{x}(\hat{\theta}_{y})\}$, indicating that the variable selection is working well.

7 Experiments with a real dataset

We show an example of auxiliary variable selection using Wine Data Set available at UCI Machine Learning Repository (Dua and Karra Taniskidou, 2017), which consists of 1 categorical variable (3 categories) and 13 continuous variables, denoted as $V_{1},\ldots,V_{13}$. For simplicity, we only use the first two categories and regard them as a latent variable $Z\in\{0,1\}$; the experiment results were similar to the other combinations. The sample size is then $n=130$ and all variables except for $Z$ are standardized. We set one of the 13 continuous variables as the observed primary variable $Y$, and set the rest of 12 variables as auxiliary variables $A_{1},\ldots,A_{12}$. For example, if $Y$ is $V_{1}$, then $A_{1},\ldots,A_{12}$ are $V_{2},\ldots,V_{13}$. The dataset is now $\{(y_{i},z_{i},a_{i,1},\ldots,a_{i,12})\}_{i=1}^{n}$, which is randomly divided into the training set with sample size $n_{tr}=86$ ($z_{i}$ is not used) and the test set with sample size $n_{te}=44$ ($a_{i,1},\ldots,a_{i,12}$ are not used).

In the experiment, we compute ${\rm AIC}_{x;b_{\ell}}$ for $B_{\ell}=(Y,A_{\ell})$, $\ell=1,\ldots,12$, and ${\rm AIC}_{x;y}$ for $Y$ from the training dataset using the model (4). We select $\hat{\theta}_{best}$ from $\hat{\theta}_{b_{1}},\ldots,\hat{\theta}_{b_{12}}$ and $\hat{\theta}_{y}$ by finding the minimum of the 13 AIC values. Thus we are selecting one of the auxiliary variables $A_{1},\ldots,A_{12}$ or not selecting any of them. It is possible to select a combination of the auxiliary variables, but we did not attempt such an experiment. For measuring the generalization error, we compute $\mathcal{L}_{x}(\hat{\theta}_{y})-\mathcal{L}_{x}(\hat{\theta}_{best})$ from the test set as

[TABLE]

where $\mathcal{D}^{te}\subset\{1,\ldots,n\}$ represents the test set. The assignment problem of $\mathcal{L}_{x}(\cdot)$ mentioned in Section 6 is avoided by a similar manner.

For each case of $Y=V_{\ell}$, $\ell=1,\ldots,13$, the above experiment was repeated $100$ times, and the experiment average of the generalization error was computed. The result is shown in Table 8. A positive value indicates that $\hat{\theta}_{best}$ performed better than $\hat{\theta}_{y}$. We observe that $\hat{\theta}_{best}$ is better than or almost the same as $\hat{\theta}_{y}$ for all cases $\ell=1,\ldots,13$, suggesting that AIC works well to select a useful auxiliary variable.

8 Conclusion

We often encounter a dataset composed of various variables. If only some of the variables are of interest, then the rest of the variables can be interpreted as auxiliary variables. Auxiliary variables may be able to improve estimation accuracy of unknown parameters but they could also be harmful. Hence, it is important to select useful auxiliary variables.

In this paper, we focused on exploiting auxiliary variables in incomplete data analysis. Usefulness of auxiliary variables is measured by a risk function based on the KL divergence for complete data. We derived an information criterion which is an asymptotically unbiased estimator of the risk function except for a constant term. Moreover, we extended a result of Stone (1977) to our setting and proved asymptotic equivalence between a variant of LOOCV and the proposed criteria. Since LOOCV requires an additional condition for its justification, the proposed criteria are preferable to LOOCV.

This study assumes that variables are different between training set and test set. There are other settings such as covariate shift (Shimodaira, 2000) and transfer learning (Pan and Yang, 2010), where distributions are different between training set and test set. It will be possible to combine these settings to construct a generalized framework. It is also possible to extend our study for taking account of a missing mechanism. We will leave these extensions as future works.

Appendix A Proofs

A.1 Proof of Lemma 1

Proof.

Taylor expansion of $\mathcal{L}_{x}(\theta)$ around $\theta=\bar{\theta}$, by formally taking it as a function of $\beta$, gives

[TABLE]

where $\nabla^{2}\mathcal{L}_{x}(\bar{\theta})=I_{x}$ is used above. By taking the expectation of both sides,

[TABLE]

where the asymptotic variance of $\hat{\beta}_{b}$ in (5) is given as

[TABLE]

∎

A.2 Proof of Lemma 2

Proof.

Taylor expansion of $\ell_{y}(\theta)$ around $\theta=\bar{\theta}$, by formally taking it as a function of $\beta$, gives

[TABLE]

where $\nabla^{2}\ell_{y}(\bar{\theta})=-I_{y}+o_{p}(1)$ is used above. By taking the expectation of both sides,

[TABLE]

In the last expression, we used (18) for the asymptotic variance of $\hat{\beta}_{b}$. For working on the second term in (19), we first derive an expression of $\hat{\beta}_{b}-\bar{\beta}$. Taylor expansion of the score function $\nabla\ell_{b}(\beta)$ around $\beta=\bar{\beta}$ gives

[TABLE]

where $\nabla^{2}\ell_{b}(\bar{\beta})=-I_{b}+o_{p}(1)$ is used above. By noticing $\nabla\ell_{b}(\hat{\beta}_{b})=0$, we thus obtain

[TABLE]

where $E[\nabla\ell_{b}(\bar{\beta})]=0$ and each term in the summation has mean zero, because $E[\nabla\log{p_{b}}(b;\bar{\beta})]=\nabla E[\log{p_{b}}(b;\bar{\beta})]=0$. Now we are back to the the second term in (19). Using (20), we have

[TABLE]

By noting $E[\nabla\ell_{b}(\bar{\beta})]=0$, the expectation of the second term in (21) is

[TABLE]

Combining (21) and (22), we have

[TABLE]

We next show that $E[\nabla^{\top}\ell_{y}(\bar{\theta})]=-\nabla^{\top}\mathcal{L}_{x}(\bar{\theta})$. Let us recall that we have assumed $q_{z|y}(z|y)=p_{z|y}(z|y;\bar{\theta})$ in (6), which leads to

[TABLE]

Therefore,

[TABLE]

Substituting this and (23) into the second term in (19), we have

[TABLE]

The first term on the right hand side in (24) is

[TABLE]

where (6) is used again in the last term. Finally (24) is rewritten as

[TABLE]

∎

A.3 Proof of Theorem 5

Proof.

First recall that we have assumed that $q_{c}(c)=p_{c}(c;\beta_{0})$, which also implies the condition (6) as $q_{z|y}(z|y)=p_{z|y}(z|y;\theta_{0})$ with $\bar{\beta}=\beta_{0}$. Thus Theorem 4 holds. Substituting $J_{b}=I_{b}$ and $K_{b,y}=I_{y}$ in the penalty term of (8), we have

[TABLE]

giving the penalty term of (10). Therefore, we only have to show (9). Noting the identity

[TABLE]

it follows from $q_{b}(b)=p_{b}(b;\beta_{0})$ that

[TABLE]

Note that the same result can be obtained from Theorem 3.3 in White (1982). Next we show $K_{b,y}=I_{y}$. Since $q_{a|y}(a|y)=p_{a|y}(a|y;\beta_{0})$,

[TABLE]

Therefore, we have

[TABLE]

∎

A.4 Proof of Theorem 6

Proof.

It follows from (14) and (15) that

[TABLE]

This and the assumption (16) imply that

[TABLE]

Under the assumption $q_{z|y}(z|y)=p_{z|y}(z|y;\bar{\theta})$,

[TABLE]

This yields that

[TABLE]

Hence, by noting $g(y;\theta)=\log p_{y}(y;\theta)+f(y;\theta)$, it holds that

[TABLE]

For evaluating the second term on the right hand side, we apply Taylor expansion to $n^{-1}\sum_{i=1}^{n}f(y_{i};\theta)$ around $\theta=\bar{\theta}$ by formally taking it as a function of $\beta$. By noting (25), this gives

[TABLE]

It follows from the law of large numbers that

[TABLE]

Hence, (18) indicates that

[TABLE]

By substituting (27) into (26), we establish that

[TABLE]

Hence, the proof completes. ∎