Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
39,001 characters · 14 sections · 23 citation commands
Uniform Inference in High-Dimensional Gaussian Graphical Models
We provide methodology and theory for uniform inference on high-dimensional graphical models with the number of target parameters being possible much larger than sample size. We demonstrate uniform asymptotic normality of the proposed estimator over $d$-dimensional rectangles and construct simultaneous confidence bands on all of the $d$ target parameters. The proposed method can be applied to test simultaneously the presence of a large set of edges in the graphical model $$X=(X_{1},\dots,X_{p})^T\sim\mathcal{N}(\mu_X,\Sigma_X).$$ Assuming that the covariance matrix $\Sigma_X$ is nonsingular, the conditional independence structure of the distribution can be conveniently represented by a graph $G = (V,E)$, where $V =\{1, \dots , p\}$ is the set of nodes and $E$ the set of edges in $V \times V$. Every pair of variables not contained in the edge set is conditionally independent given all remaining variables. If the vector $X$ is normally distributed, every edge corresponds to a non-zero entry in the inverse covariance matrix (Lauritzen (1996)) lauritzen1996graphical.\\ \\ In the last decade, significant progress has been made on estimation of a large precision matrix in order to analyze the dependence structure of a high-dimensional normal distributed random variable. There are mainly two common approaches to estimate the entries of a precision matrix. The first approach is a penalized likelihood estimation approach with a lasso-type penalty on entries of the precision matrix, typically referred to as the graphical lasso. This approach has been studied in several papers, see e.g Lam and Fan (2009) lam2009sparsistency, Rothman et al. (2008) rothman2008sparse, Ravikumar et al. (2011) ravikumar2011high and Yuan and Lin (2007) yuan2007model. The second approach, first introduced by Meinshausen and B\"uhlmann (2006) meinshausen2006high, is neighborhood based. It estimates the conditional independence restrictions separately for each node in the graph and is hence equivalent to variable selection for Gaussian linear models. The idea of estimating the precision matrix column by column by running a regression for each variable against the rest of variables was further studied in Yuan (2010) yuan2010high, Cai, Liu and Zhou (2011) cai2011constrained and Sun and Zhang (2013) sun2013sparse. \\ \\ In this paper, we do not aim to estimate the whole precision matrix but we focus on quantifying the uncertainty of recovering its support by providing a significance test for a set of potential edges. In recent years, statistical inference for the precision matrix in high-dimensional settings has been studied, e.g in Jankov\'{a} and van de Geer (2016) jankova2017honest and Ren et al. (2015) ren2015asymptotic. Both approaches lead to an estimate that is elementwise asymptotically normal and enables testing for low-dimensional parameters of the precision matrix using standard procedures such as Bonferroni-Holm correction.\\ In contrast to these existing results, our method explicitly allows for testing a joint hypothesis without correction for multiple testing and conducting inference for a growing number of parameters using high dimensional central limit results. In particular, our results rely on approximate sparsity instead of row sparsity which restricts the number of non-zero entries of each row of the precision matrix to be at most $s\ll n$ that is in many applications a questionable assumption. In order to provide theoretical results, fitting the problem of support discovery in Gaussian graphical models into a general Z-estimation setting with a high-dimensional nuisance function is key. Inference on a (multivariate) target parameter in general Z-estimation problems in high dimensions is covered in Belloni et al. (2014) belloni2014uniform, Belloni et al. (2018) belloni2018uniformly and Chernozhukov et al. (2017) chernozhukov2017double. To conduct inference on a high-dimensional target parameter, uniform estimation rates and sparsity guarantees of the nuisance function are crucial. In this context, we formally apply recent results from Belloni et al. (2018) belloni2018uniformly to ensure sufficient fast convergence rate of the lasso estimator under approximate sparsity conditions. Moreover, we provide auxiliary results for the square-lasso estimator establishing uniform estimation rates and sparsity guarantees in a random design under approximate sparsity conditions that might be of independent interest for related problems in high-dimensional linear models.
The rest of this paper is organized as follows. In Section (ref), we formally define the setting and introduce the notation that will be used fitting the problem of support discovery in Gaussian graphical models into a general Z-estimation problem with a high-dimensional nuisance function. In Section (ref), we outline the estimation procedure of the high-dimensional target parameter and the conditions that are needed to achieve our main theorem presented in Section (ref). Section (ref) provides implementation details and shows how our estimation procedure can be modified by cross-fitting to improve small sample properties. Section (ref) provides a simulation study on the proposed method. The supplementary material includes additional technical material. The proof of our main theorem is provided in Appendix (ref). The uniform nuisance function estimation is discussed in Appendix (ref). Appendix (ref) formally discusses conditions for the uniform convergence rates of the lasso estimator. Finally, Appendix (ref) provides auxiliary results for the square-lasso estimator.
Let $$X=(X_{1},\dots,X_{p})^T\sim\mathcal{N}(\mu_X,\Sigma_X)$$ be a $p$-dimensional random variable. For all $(j,k)\in E$ with $j\neq k$, assume that $$X_j=\sum\limits_{\substack{l =1\\ l\neq j}}^p \beta_{l}^{(j)}X_l+\varepsilon^{(j)}=\beta^{(j)} X_{-j}+\varepsilon^{(j)}$$ and $$X_k=\gamma^{(j,k)}X_{-\{j,k\}}+\nu^{(j,k)},$$ where $\mathbb{E} [\varepsilon^{(j)}|X_{-j}]=0$ and $\mathbb{E} [X_{-\{j,k\}}\nu^{(j,k)}]=0$. Define the column vector $$\Gamma^{(j)}=\left(-\beta^{(j)}_1,\dots,-\beta^{(j)}_{j-1},1,-\beta^{(j)}_{j+1},\dots,-\beta^{(j)}_{p}\right)^T.$$ One may show
where $\Phi_0^{j}$ is the $j$-th column of the precision matrix $\Phi_0=\Sigma_X^{-1}$ jankova2017honest. Hence
for all $j\neq k$. Assume that we are interested in the following set of potential edges $$\mathcal{M}:=\{m_1,\dots,m_{d_n}\}$$ where the number of edges $d_n$ may increase with sample size $n$. In the following the dependence on $n$ is omitted to simplify the notation. In order to test whether all variables $X_j$ and $X_k$ are conditionally independent with $m_r=(j_r,k_r)$ for a $r\in\{1,\dots,d\}$, we have to estimate our target parameter $$\theta_0=(\theta_{m_{1}},\dots,\theta_{m_{d}})^T:=(\beta^{(j_1)}_{k_1},\dots,\beta^{(j_d)}_{k_d})^T.$$ The setting above fits in the general Z-estimation problem of the form $$\mathbb{E} \left[\psi_{m_r}\big(X,\theta_{m_r},\eta_{m_r}\big)\right]=0$$ for all $r=1,\dots,d$ with nuisance parameters $$\eta_{m_r}=\left(\beta^{(j)}_{-k},\gamma^{(j,k)}\right)$$ where $\beta^{(j)}_{-k}\equiv\beta^{(m_r)}$ and $\gamma^{(j,k)}\equiv\gamma^{(m_r)}$. The score functions are defined by
for $m_r=(j_r,k_r)\equiv (j,k)$, $\eta=(\eta^{(1)},\eta^{(2)})$ and $r=1,\dots,d$. Without loss of generality we assume $j>k$ for all tuples $m_r\in \mathcal{M}$.
Let $X^{(i)}$, $i=1,\dots,n$ be i.i.d. random vectors.\\ At first we estimate the nuisance parameter $\eta_{m_r}=\big(\eta_{m_r}^{(1)},\eta_{m_r}^{(2)}\big)$ by running a lasso/ post-lasso/ square-root lasso regression of $X_j$ on $X_{-j}$ to compute $(\tilde{\theta}_{m_r},\hat{\eta}_{m_r}^{(1)})$ and a lasso/ post-lasso/ square-root lasso regression of $X_k$ on $X_{-m_r}$ to compute $\hat{\eta}_{m_r}^{(2)}$ for each $(j,k)=m_r\in\mathcal{M}$. The estimator $\hat{\theta}_0$ of the target parameter $$\theta_0=(\theta_{m_1},\dots,\theta_{m_{d}})^T$$ is defined as the solution of
where $\epsilon_{n}=o\left(\delta_nn^{-1/2}\right)$ is the numerical tolerance and $(\delta_n)_{n\ge 1}$ a sequence of positive constants converging to zero.\\ \\ Assumptions A1-A4.\\ Let $a_n:=\max(d,p,n,e)$ and $C$ a strictly positive constant independent of $n$ and $r$. The following assumptions hold uniformly in $n\ge n_0,P\in\mathcal{P}_n$:
The condition (ref) is a standard approximate sparsity condition that is discussed in more detail in comment (ref). The number of relevant variables $s_n\equiv s$ captured by the regression coefficient $\beta^{(1,m_r)}$ respectively $\gamma^{(1,m_r)}$ can grow with the sample size. The coefficient $\beta^{(2,m_r)}$ respectively $\gamma^{(2,m_r)}$ is the approximate sparse part of the true regression coefficient. This misspecification of a sparse model is controlled by condition (ref). The growth condition (ref) ensures that $s^2\log^4(a_n)/n$ converges towards zero with at least polynomial speed. If this convergence is too slow ($\tilde{q}\ge 9$) the condition on the growth rate of the number of tested edges become more restrictive. In general, both the number of parameters $p$ and the number of relevant variables $s$ can grow with the sample size in a balanced way. If $s$ is fixed, the number of potential parameters $p$ can grow at an exponential rate with the sample size. This means that the set of potential variables can be much larger than the sample size, only the number of relevant variables $s$ has to be smaller than the sample size. This situation is common for Lasso-based estimators. Condition (ref) restricts the parameter spaces and ensures that the true coefficients are well behaved. The condition (ref) is a standard eigenvalue condition that restricts the correlation between the components of $X$ and bounds the variances of each $X_j$ from below and above. Assumptions (ref)-(ref) combined with the normal distribution of $X$ imply the conditions (ref)-(ref) from theorem (ref) which enables us to estimate the nuisance parameter sufficiently fast by lasso and post-lasso. To ensure a sufficiently fast convergence rate and sparsity guarantees of the square-root lasso estimator further model assumptions are needed.
We will prove that the assumptions of Corollary $2.2$ from Belloni et al. (2018) belloni2018uniformly hold and hence we are able to use their results to construct confidence intervals even for a growing number of hypothesis $d=d_n$. Define
and the corresponding estimators
for $r=1,\dots,d$. To construct confidence intervals we will employ the Gaussian multiplier bootstrap. Define $$\hat{\psi}_{m_r}(X):=-\hat{\sigma}_{m_r}^{-1}\hat{J}_{m_r}^{-1}\psi_{m_r}(X,\hat{\theta}_{m_r},\hat{\eta}_{m_r})$$ and the process $$\hat{\mathcal{N}}:=\left(\hat{\mathcal{N}}_{m_r}\right)_{m_r\in\mathcal{M}}=\left(\frac{1}{\sqrt{n}}\sum\limits_{i=1}^n\xi_i\hat{\psi}_{m_r}\big(X^{(i)}\big)\right)_{m_r\in\mathcal{M}}$$ where $(\xi_i)_{i=1}^n$ are independent standard normal random variables which are independent from $\big(X^{(i)}\big)_{i=1}^n$. We define $c_{\alpha}$ as the $(1-\alpha)$-conditional quantile of $\sup_{m_r\in\mathcal{M}}|\hat{\mathcal{N}}_{m_r}|$ given the observations $\big(X^{(i)}\big)_{i=1}^n$. The following theorem is the main result of our paper and establishes simultaneous confidence bands for the target parameter $\theta_0$.
\ \\ Using theorem (ref) we are able to construct standard confidence regions which are uniformly valid over a large set of variables and we can check null hypothesis of the form: $$H_0: \mathcal{M}\cap E = \emptyset.$$
We implemented a function that will be added to the $R$-package $hdm$ and estimates the target coefficients $$(\theta_{m_{1}},\dots,\theta_{m_{d}})^T=(\beta^{(j_1)}_{k_1},\dots,\beta^{(j_d)}_{k_d})^T$$ corresponding the considered set of potential edges $$\mathcal{M}:=\{m_1,\dots,m_{d_n}\}$$ by the proposed method described in section (ref). It can be used to perform hypothesis tests with asymptotic level $\alpha$ based on the different confidence regions described in comment (ref). The nuisance function can be estimated by lasso, post-lasso or square-root lasso.
In general Z- estimation problems where a so called debiased or double machine learning (DML) method is used to construct confidence intervals, it is common to use cross-fitting in order to improve small sample properties. A detailed discussion of cross-fitted DML can be found in Chernozhukov et al. (2017) chernozhukov2017double. The following algorithm generalizes our proposed method to a $K$-fold cross fitted version. We assume that $n$ is divisible by $K$ in order to simplify notation.
The confidence region above corresponds to ((ref)). Confidence regions corresponding to ((ref)) or ((ref)) can be constructed in an analogous way.
This section provides a simulation study on the proposed method. In each example the precision matrix of the Gaussian graphical model is generated as in the $R$-package $huge$ zhao2012huge. Hence, the corresponding adjacency matrix $A$ is generated by setting the nonzero off-diagonal elements to be one and each other element to be zero. To obtain a positive definite pre-version of the precision matrix we set $$\Phi_{pre}:= v\cdot A+(|\Lambda_{\min}(v\cdot A)|+0.1+u)\cdot I_{p\times p}.$$ Here $v=0.3$ and $u=0.1$ are chosen to control the magnitude of partial correlations. The covariance matrix $\Sigma$ is generated by inverting $\Phi_{pre}$ and scaling the variances to one. The corresponding precision matrix $\Phi$ is given by $\Sigma^{-1}$. For a given $p$ we generate $n=200$ independent samples of $$X=(X_1,\dots,X_p)\sim\mathcal{N}(0,\Sigma)$$ and evaluate whether our test statistic would reject the null hypothesis for a specific set of edges $\mathcal{M}$ which satisfies the null hypothesis. Finally the acceptance rate is calculated over $l=1000$ independent simulations for a given confidence level $1-\alpha=0.95$.
In our simulation study we estimate the correlation structure of four different designs that are described in the following.
Each pair of off-diagonal elements of the covariance matrix of the first $p-1$ regressors is randomly set to non-zero with probability $prob = 5/p$. The last regressor is added as an independent random variable. It results in about $(p-1)\cdot(p-2)\cdot prob /2$ edges in the graph. The corresponding precision matrix is of the form $$ \Phi:=\left(
\right) $$ where $B$ is a sparse matrix. We test the hypothesis, whether the last regressor is independent from all other regressors, corresponding to $$\mathcal{M}=\{(p,1),\dots,(p,p-1)\}.$$
The regressors are evenly partitioned into $g=4$ disjoint groups. Each pair of off-diagonal elements $\Phi_{(i,j)}$ is set non-zero with probability $prob=5/p$, if both $i$ and $j$ belong to the same group. It results in about $g\cdot(p/g)\cdot(p/g-1)\cdot prob/2$ edges in the graph. The precision Matrix is of the form $$ \Phi:=\left(
\right) $$ where each block $B_i$ is a sparse matrix. We test the hypothesis that the first two hubs are conditionally independent. This corresponds to testing the tuples $$\mathcal{M}=\{(1,p/4+1),\dots,(1,p/2),(2,p/4+1),\dots,(p/4,p/2)\}.$$
In this example we generate a random graph structure as in example $1$, but instead of setting the other elements of the adjacency matriy $A$ to zero we generate independent random entries from a uniform distribution on $[-a,a]$ with $a=1/20$. This results in a precision matrix of the form $$ \Phi:=\left(
\right) $$ where $B$ is not a sparse matrix anymore. We then again test the hypothesis, whether the last regressor is independent from all other regressors, corresponding to $$\mathcal{M}=\{(p,1),\dots,(p,p-1)\}.$$
By setting $$\Phi:=I_{p\times p}$$ we generate samples of $p$ independent normal distributed random variables. We can test the hypothesis whether the regressors are independent by choosing $$\mathcal{M}=\{(1,2),\dots,(1,p),(2,3),\dots,(p-1,p)\}.$$
We provide simulated acceptance rates of our proposed estimation procedure with $B=1000$ bootstrap samples for all of the examples above. Confidence Intervall I corresponds to the standard case in ((ref)), whereas Confidence Intervall II is based on the approximation of the sphere in ((ref)). In summary, the results reveal that the empirical acceptance rate is, on average, close to the nominal level of $95\%$ with a mean absolute deviation of $2.581\%$ over all simulations. The Confidence Intervall II has got a mean absolute deviation of $1.875\%$ and performs significantly better than Confidence Intervall I with a mean absolute deviation of $3.287\%$. More complex S-sparsely convex sets seem to result in better acceptance rates, whereas higher exponents do not improve the rates. The lowest mean absolute deviation ($1.138\%$) is achieved in table 2 for $S=5$, $exp=1$ and without cross-fitting.