Sensitivity analysis for inference with partially identifiable covariance matrices

Abstract

In some multivariate problems with missing data, pairs of variables exist that are never observed together. For example, some modern biological tools can produce data of this form. As a result of this structure, the covariance matrix is only partially identifiable, and point estimation requires that identifying assumptions be made. These assumptions can introduce an unknown and potentially large bias into the inference. This paper presents a method based on semidefinite programming for automatically quantifying this potential bias by computing the range of possible equal-likelihood inferred values for convex functions of the covariance matrix. We focus on the bias of missing value imputation via conditional expectation and show that our method can give an accurate assessment of the true error in cases where estimates based on sampling uncertainty alone are overly optimistic.

Appendix: Algorithmic details

From Sect. 3.2, the inner optimization problem we wish to solve is

$$\begin{aligned}&\underset{\varSigma }{\mathrm{minimize}} \;\; -\frac{1}{t}\log |\varSigma | + \varSigma _{j,J_k}\hat{\varSigma }_{J_k,J_k}^{-1}(x_{J_k}-\hat{\mu }_{J_k})\\&\text {subject to } \;\; \varSigma _{ab} = \hat{\varSigma }_{ab} \text { for } (a,b)\in \bigcup _{\ell } \left( J_{\ell }\times J_{\ell }\right) . \end{aligned}$$

This optimization problem, with fixed $t$, can be solved by generalized gradient descent. The gradient of the objective is

$$\begin{aligned} \nabla (\mathrm{objective})&= -\frac{1}{t} \varSigma ^{-1} + C \\ C_{ab}&= \left\{ \begin{array}{ll} \left( \hat{\varSigma }_{J_k,J_k}^{-1}(x_{J_k}-\hat{\mu }_{J_k})\right) [b] &{} \quad a=j \text { and } b\in J_k\\ \left( \hat{\varSigma }_{J_k,J_k}^{-1}(x_{J_k}-\hat{\mu }_{J_k})\right) [a] &{} \quad a\in J_k \text { and } b=j\\ 0 &{} \quad \text {otherwise} \end{array}\right. \end{aligned}$$

If we initialize the first time with $\varSigma = \hat{\varSigma } \in \mathcal S$, we want to remain within $\mathcal S$ with each step. Therefore, we project the gradient into the linear space $\{\varSigma : \varSigma _{ij} = \hat{\varSigma }_{ij} \forall (i,j) \in \bigcup _k \left( J_k\times J_k\right) \}$. This is equivalent to holding the coordinates $(i,j) \in \bigcup _k \left( J_k\times J_k\right)$ fixed, and only taking gradient steps in the other directions.

With a step size of $\delta$, the update becomes

$$\begin{aligned} \begin{array}{ll} \varSigma ^{(t+1)}_{ab} = \varSigma ^{(t)}_{ab} + \delta \left( \frac{1}{t} (\varSigma ^{(t)})_{ab}^{-1} - C^{(t)}_{ab}\right) &{} \quad (a,b) \notin \bigcup _{\ell } \left( J_{\ell }\times J_{\ell }\right) \\ \varSigma ^{(t+1)}_{ab} = \varSigma ^{(t)}_{ab} &{} \quad (a,b) \in \bigcup _{\ell } \left( J_{\ell }\times J_{\ell }\right) \end{array} \end{aligned}$$

Using warm starts, we repeatedly solve this problem with increasing $t$ to obtain the final solution. This barrier method is discussed in Boyd and Vandenberghe (2004), which recommends using a sequence of $t$ that increase by a factor of $\mu$ (around 10–20) at each outer loop iteration. More details of the method can be found in Section 11.3.1 of that book.

We also include acceleration, as in Banerjee et al. (2006) and Beck and Teboulle (2010), among others. This is shown in Fig. 5 to give practically significant improvements in algorithm timings. The final algorithm is shown in Algorithm 1.

Some care needs to be taken in selecting the step size $\delta$. We choose it by backtracking to result in a decrease in the objective, to remain inside the positive semidefinite cone, and to satisfy the majorization requirements of generalized gradient descent (Beck and Teboulle 2010). This sub-algorithm is shown as Algorithm 2.