\section{Introduction}
Consider a regression problem with $d$ independent covariates $\mathsf{x}_1,\ldots,\mathsf{x}_d$ and a binary outcome variable $\mathsf{y}$. The covariates are assumed to be uniformly distributed over $\{+1,-1\}$, and the conditional probabilities of the outcome given the covariates are assumed to be logistic:
\begin{align}
\mathrm{Pr}\left(\mathsf{y}=+1| \mathsf{x}_V=x_V\right) & =\sigma\Big(\sum_{\{i,j\}\in E}\beta_{\{i,j\}}x_ix_j\Big),\label{conditional1}\\
\mathrm{Pr}\left(\mathsf{y}=-1| \mathsf{x}_V=x_V\right) & =\sigma\Big(-\sum_{\{i,j\}\in E}\beta_{\{i,j\}}x_ix_j\Big)\label{conditional2}
\end{align}
for some real constants $\beta_{\{i,j\}}$'s, where $V:=\{1,\ldots,d\}$, $E:=\{\{i,j\}:i,j\in V, i \neq j\}$, $\mathsf{x}_V:=(\mathsf{x}_i:i\in V)$, and $\sigma(x):=e^{x}/(1+e^x)$ is the sigmoid function. It is straightforward to verify that $\sigma(x)+\sigma(-x)=1$ for any $x \in \mathbb{R}$, so we have $\mathrm{Pr}\left(\mathsf{y}=+1| \mathsf{x}_V=x_V\right)+\mathrm{Pr}\left(\mathsf{y}=-1| \mathsf{x}_V=x_V\right)=1$ for any $x_V\in \{+1,-1\}^d$.

For any two distinct $i,j\in V$, we say that the covariates $\mathsf{x}_i$ and $\mathsf{x}_j$ {\em interact} if and only if $\beta_{\{i,j\}} \neq 0$. Let $G=(V,I)$ be a simple graph with the vertex set $V$ and edge set $I:=\{\{i,j\} \in E: \beta_{\{i,j\}}\neq 0\}$. Then $G$ captures all pairwise interactions between the covariates in determining the odds of the outcome of interest. Our goal is to recover the graph $G$ from independently identically distributed (i.i.d.) samples of $(\mathsf{x}_V,\mathsf{y})$.

Our main motivation for considering the above pairwise interaction problem is from computational biology, where each covariate $\mathsf{x}_i$ represents the expression of a biomarker (a gene or an environmental factor), and the variable $\mathsf{y}$ represents the phenotypic outcome with respect to a specific phenotype. In computational biology, many complex diseases, such as cancer and diabetes, are conjectured to have complicated underlying disease mechanisms~\cite{mo2003,an2007,ed2010,ch2007,ho2011,adl2013,sa2015}. Multiple candidate risk factors, either genetic or environmental, and their interactions are known to play critical roles in triggering and determining the development of a large family of diseases~\cite{mo2003,an2007,ed2010,ch2007,ho2011,adl2013,sa2015}. Identifying interactive effects among profiled variables not only helps with more accurate identification of critical risk factors for outcome prediction, but also helps reveal functional interactions and understand aberrant system changes that are specifically related to the outcome for effective systems intervention. Our model, of course, is a simplified one from real-world situations, and is studied here since it captures some essential features of the problem (as we shall see shortly) while being relatively simple.

Note that if we let
\begin{align}
\mathsf{z}_{\{i,j\}}:=\mathsf{x}_i\mathsf{x}_j, \quad \forall \{i,j\}\in E
\end{align}
and consider $\mathsf{z}_{\{i,j\}}$'s (instead of $\mathsf{x}_i$'s) as the covariates, then the problem of recovering the graph $G$ can be viewed as a {\em feature selection} problem in statistics and machine learning. In \cite{sa2007}, a basic approach for feature selection is to first use Shannon's mutual information \cite{co2006} to measure the {\em marginal} effects of the covariates on the outcome, and then select the features based on the ranking of the mutual information. More advanced approaches such as the immensely popular mRMR method \cite{pe2005} make incremental selections while taking into account both the relevance to the outcome and the redundancy among the selected features. However, even though Shannon's mutual information provides a compact and model-free measure of correlation between the covariates and the outcome, which is well accepted in the statistics and computational biology communities, it is a complex function of the underlying joint distribution and hence difficult to analyze and estimate from limited data samples. As a result, when applied to specific regression models, the performance of the generic feature selection algorithms is usually difficult to characterize.

Motivated by the recent progress on learning Ising models over arbitrary graphs \cite{br2005}, in this paper we propose a quantity called {\em influence} as a measure of the marginal effects of $\mathsf{z}_{\{i,j\}}$'s on the outcome $\mathsf{y}$. Compared with Shannon's mutual information, influence is a simple function of the low-order joint probabilities between $\mathsf{x}_i$'s and the outcome $\mathsf{y}$, and hence is much easier to analyze and estimate. When the underlying graph $G$ is known to be {\em acyclic}, we show that, a simple algorithm that is based on a maximum-weight spanning tree with respect to the ``plug-in" estimate of the influences and followed by simple thresholding operations, not only has strong theoretical performance guarantees, but can also outperform the generic feature selection algorithms for recovering $G$ from i.i.d. samples of $(\mathsf{x}_V,\mathsf{y})$.

The rest of the paper is organized as follows. In Section~\ref{sec:Id}, we show that any acyclic $G$ can be identified from the influences of $\mathsf{z}_{\{i,j\}}$'s on the outcome $\mathsf{y}$. Building on the results from Section~\ref{sec:Id}, in Section~\ref{sec:Det} we show that any acyclic $G$ can be recovered with probability at least $1-\epsilon$ from $n$ i.i.d. samples of $(\mathsf{x}_V,\mathsf{y})$, where $n=\Theta\left(d\log(d^2/\epsilon)\right)$. In Section~\ref{sec:Ext}, we extend our results of the above sections to the model involving both individual effects and cooperative interactions. In Section~\ref{sec:Sim}, we use computer simulations to demonstrate that the proposed algorithm can outperform the generic feature selection algorithms. Finally, in Section~\ref{sec:Con} we conclude the paper with some remarks.

{\em Notation}. Random variables are written in serif font, and sets are written in capital letters.

\section{Identification of Cooperative Interactions from Low-Order Joint Probabilities}\label{sec:Id}
Our main result in this section is to show that any {\em acyclic} $G$ can be identified from the low-order joint probabilities $(p(\mathsf{x}_i,\mathsf{x}_j,\mathsf{y}):\{i,j\}\in E)$. Towards this goal, let $w$ be a weight assignment over $E$ given by:
\begin{align}
w_{\{i,j\}} :&= \left|\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1\right)-\right.\nonumber\\
& \hspace{40pt} \left.\mathrm{Pr}\left(\mathsf{y}=-1|\mathsf{x}_i=+1,\mathsf{x}_j=+1\right)\right|\label{eq:W}\\
& = \left|2\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1\right)-1\right|\label{eq:W2}\\
&=\left|8\mathrm{Pr}\left(\mathsf{x}_i=+1,\mathsf{x}_j=+1,\mathsf{y}=+1\right)-1\right|\label{eq:W3}
\end{align}
for any $\{i,j\}\in E$. Here, \eqref{eq:W2} follows from the fact that
\begin{align*}
\mathrm{Pr}&\left(\mathsf{y}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1\right)+\\
&\mathrm{Pr}\left(\mathsf{y}=-1|\mathsf{x}_i=+1,\mathsf{x}_j=+1\right)=1
\end{align*}
and \eqref{eq:W3} is due to the fact that $\mathrm{Pr}(\mathsf{x}_i=+1,\mathsf{x}_j=+1)=1/4$.

The following proposition helps to clarify the meaning of the weight assignment $w$ as defined in \eqref{eq:W}.

\begin{prop}[Influence]
\label{prop:Inf}
Assume that $G=(V,I)$ is acyclic. We have
\begin{align}
w_{\{i,j\}} &= \left|\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1\right)-\right.\nonumber\\
& \hspace{40pt} \left.\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=-1\right)\right|\label{eq:Inf}
\end{align}
for any $\{i,j\}\in E$.
\end{prop}

\begin{proof}
See Section~\ref{pf:Inf}.
\end{proof}

By \eqref{eq:Inf}, $w_{\{i,j\}}$ indicates whether the product $\mathsf{z}_{\{i,j\}}=\mathsf{x}_i\mathsf{x}_j$ has any {\em influence} on the event $\mathsf{y}=+1$ and hence can be a useful indication on whether $\{i,j\} \in I$. This intuition is partially justified by the following proposition.

\begin{prop}[Direct influence]
\label{prop:DI}
Assume that $G=(V,I)$ is acyclic. We have $w_{\{i,j\}}>0$ for any $\{i,j\} \in I$.
\end{prop}

\begin{proof}
See Section~\ref{pf:DI}.
\end{proof}

We say that the product $\mathsf{z}_{\{i,j\}}=\mathsf{x}_i\mathsf{x}_j$ has a {\em direct influence} on the outcome $\mathsf{y}$ if $\{i,j\}\in I$. The above proposition guarantees that direct influences are strictly positive when $G$ is acyclic. The following proposition provides a {\em partial} converse to Proposition~\ref{prop:DI}.

\begin{prop}[Zero influence]
\label{prop:ZI}
Assume that $G=(V,I)$ is acyclic. Then for any two distinct $i,j\in V$, we have $w_{\{i,j\}}=0$ if $i$ and $j$ are disconnected in $G$, or the unique path between $i$ and $j$ in $G$ has an even length.
\end{prop}

\begin{proof}
See Section~\ref{pf:ZI}.
\end{proof}

\begin{thm}[Union of stars]
Assume that each connected component of $G=(V,I)$ is a star. Then for any two distinct $i,j \in V$, we have $\{i,j\} \in I$ if and only if $w_{\{i,j\}}>0$.
\end{thm}

\begin{proof}
This follows immediately from Propositions~\ref{prop:DI} and \ref{prop:ZI} and the fact that if each connected component of $G$ is a star (which implies that $G$ is acyclic), then any two distinct $i,j\in V$ such that $\{i,j\}\not\in I$ must be either disconnected (if they belong to two different connected components) or connected by a unique path of length two (if they belong to the same connected component) in $G$.
\end{proof}

The following example, however, shows that the converse of Proposition~\ref{prop:DI} is {\em not} true in general. Consider $d=4$, $I=\{\{1,2\},\{2,3\},\{3,4\}\}$, and $\beta_{\{1,2\}}=\beta_{\{2,3\}}=\beta_{\{3,4\}}=1$. Note that the graph $G=(V,I)$ here is acyclic and the unique path between $1$ and $4$ is of length three. It is straightforward to calculate that
\begin{align*}
w_{\{1,4\}}&=\frac{(e-1)^3}{2(e^3+1)}>0,
\end{align*}
even though $\{1,4\} \not\in I$.

For any $\{i,j\}\in E\setminus I$, we say that the product $\mathsf{z}_{\{i,j\}}=\mathsf{x}_i\mathsf{x}_j$ has an {\em indirect influence} on the outcome $\mathsf{y}$ if $w_{\{i,j\}}>0$. Due to the possible existence of indirect influences, unlike the unions of stars, a general acyclic $G$ cannot be recovered via edge-by-edge identifications.

The following proposition, however, shows that indirect influences are {\em locally} dominated by direct influences.

\begin{prop}[Indirect influence]
\label{prop:II}
Assume that $G=(V,I)$ is acyclic, and let $\left\{\{i_1,i_2\},\{i_2,i_3\},\ldots,\{i_m,i_{m+1}\}\right\}$ be a path of length $m \geq 2$ in $G$. Then, we have $w_{\{i_1,i_{m+1}\}}<w_{\{i_s,i_{s+1}\}}$ for any $s\in\{1,\ldots,m\}$.
\end{prop}

\begin{proof}
Note that when $m$ is even, by Propositions~\ref{prop:DI} and \ref{prop:ZI} we have $w_{\{i_1,i_{m+1}\}}=0<w_{\{i_s,i_{s+1}\}}$ for any $s\in\{1,\ldots,m\}$. Therefore, we only need to consider the cases where $m$ is odd, for which the proof can be found in Appendix~\ref{pf:II}.
\end{proof}

Let $D:=\{\{i,j\}\in E: \mbox{$i$ and $j$ are disconnected in $G$}\}$. A weight assignment $u$ over $E$ is said to have {\em strict separation} between $I$ and $D$ if there exists a real constant $\eta\geq0$ such that $u_{\{i,j\}} > \eta \geq u_{\{k,l\}}$ for any $\{i,j\}\in I$ and $\{k,l\}\in D$. The consequence of strict separation and local dominance is summarized in the following proposition.

\begin{prop}\label{prop:Cond}
Assume that $G=(V,I)$ is acyclic, let $u$ be a weight assignment over $E$ satisfying: 1) (strict separation) there exists a real constant $\eta\geq0$ such that $u_{\{i,j\}} > \eta \geq u_{\{k,l\}}$ for any $\{i,j\}\in I$ and $\{k,l\}\in D$; and 2) (local dominance) $u_{\{i_1,i_{m+1}\}}<u_{\{i_s,i_{s+1}\}}$ for any path $\left\{\{i_1,i_2\},\{i_2,i_3\},\ldots,\{i_m,i_{m+1}\}\right\}$ of length $m \geq 2$ in $G$ and any $s\in\{1,\ldots,m\}$. Then for any maximum-weight spanning tree $G'=(V,T)$ with respect to the weight assignment $u$, we have $I=T\cap W$ where $W:=\{\{i,j\}\in E: u_{\{i,j\}} > \eta\}$.
\end{prop}

\begin{proof}
See Section~\ref{pf:Cond}.
\end{proof}

The following theorem is the main result of this section.

\begin{thm}[Acyclic graphs]\label{thm:AG}
Assume that $G=(V,I)$ is acyclic, and let $G'=(V,T)$ be a maximum-weight spanning tree with respect to the weight assignment $w$ as defined in \eqref{eq:W}. Then, for any two distinct $i,j\in V$ we have $\{i,j\}\in I$ if and only if $\{i,j\}\in T$ and $w_{i,j}>0$.
\end{thm}

\begin{proof}
Note that by Propositions~\ref{prop:DI} and \ref{prop:ZI}, we have $w_{\{i,j\}} > 0=w_{\{k,l\}}$ for any $\{i,j\}\in I$ and $\{k,l\}\in D$. By Proposition~\ref{prop:II}, we have $w_{\{i_1,i_{m+1}\}}<w_{\{i_s,i_{s+1}\}}$ for any path $\left\{\{i_1,i_2\},\{i_2,i_3\},\ldots,\{i_m,i_{m+1}\}\right\}$ of length $m \geq 2$ in $G$ and any $s\in\{1,\ldots,m\}$. The theorem thus follows directly from Proposition~\ref{prop:Cond} with $u=w$ and $\eta=0$.
\end{proof}

\section{Detection of Cooperative Interactions from Finite Data Samples}\label{sec:Det}
Let $(\mathsf{x}_V[t],\mathsf{y}[t])$, $t=1,\ldots,n$ be $n$ i.i.d. samples of $(\mathsf{x}_V,\mathsf{y})$. To recover the graph $G$, we shall assign to each $\{i,j\} \in E$ a weight that is based on the {\em empirical} joint probability:
\begin{align}\label{eq:What}
\hat{\mathsf{w}}_{\{i,j\}}&:=\left|\frac{8}{n}\sum_{t=1}^n\mathbbm{1}((\mathsf{x}_i[t],\mathsf{x}_j[t],\mathsf{y}[t])=(+1,+1,+1))-1\right|,
\end{align}
where $\mathbbm{1}(\cdot)$ is the indicator function. By \eqref{eq:W3}, for any $\{i,j\} \in E$ we have $\hat{\mathsf{w}}_{\{i,j\}}$ converging to $w_{\{i,j\}}$ in probability in the limit as $n \rightarrow \infty$. The following simple proposition, which follows directly from the well-known Hoeffding's inequality \cite{ho1963}, provides a bound on the rate at which the weight assignment $\hat{\mathsf{w}}$ converges {\em uniformly} to $w$.

\begin{prop}\label{prop:Chernoff}
For any $\{i,j\}\in E$ and $\eta>0$,
\begin{align}
&\mathrm{Pr}\Big(\left|\hat{\mathsf{w}}_{\{i,j\}}-w_{\{i,j\}}\right|\ge \eta\Big) \leq 2e^{-n\eta^2/32}.
\end{align}
\end{prop}

\begin{proof}
Note that \begin{align*}
\Big|\hat{\mathsf{w}}_{\{i,j\}}&-w_{\{i,j\}}\Big|
\le 8\bigg|\mathrm{Pr}\left(\mathsf{x}_i=+1,\mathsf{x}_j=+1,\mathsf{y}=+1\right)\\
&-\frac{1}{n}\sum_{t=1}^n\mathbbm{1}((\mathsf{x}_i[t],\mathsf{x}_j[t],\mathsf{y}[t])=(+1,+1,+1))\bigg|.
\end{align*}
We then finish the proof by applying the Hoeffding's inequality \cite{ho1963}.
\end{proof}


The following propositions, which are generalizations of Propositions~\ref{prop:DI} and \ref{prop:II} respectively, play a key role in adapting the results of Theorem~\ref{thm:AG} from the weight assignment $w$ to $\hat{\mathsf{w}}$.

\begin{prop}\label{prop:PC1}
Assume that $G=(V,I)$ is acyclic and for any $\{i,j\}\in I$ we have $|\beta_{\{i,j\}}|\in [\lambda,\mu]$ for some $\mu \geq \lambda >0$. Let
\begin{align}
\gamma:=\sqrt{\frac{2}{\pi d}} \left[\sigma(\lambda+3\mu)-\sigma(-\lambda+3\mu)\right]>0.\label{eq:gamma}
\end{align}
We have $w_{\{i,j\}}\geq\gamma$ for any $\{i,j\} \in I$.
\end{prop}

\begin{proof}
See Section~\ref{pf:DI}.
\end{proof}

\begin{prop}\label{prop:PC2}
Assume that $G=(V,I)$ is acyclic and for any $\{i,j\}\in I$ we have $|\beta_{\{i,j\}}|\in [\lambda,\mu]$ for some $\mu \geq \lambda >0$. We have $w_{\{i_1,i_{m+1}\}}\leq w_{\{i_s,i_{s+1}\}}-\gamma$ for any path $\left\{\{i_1,i_2\},\{i_2,i_3\},\ldots,\{i_m,i_{m+1}\}\right\}$ of length $m \geq 2$ in $G$ and any $s\in\{1,\ldots,m\}$, where $\gamma$ is defined as in \eqref{eq:gamma}.
\end{prop}

\begin{proof}
See Section~\ref{pf:II}.
\end{proof}

Given Propositions~\ref{prop:PC1} and \ref{prop:PC2}, it is clear that if the estimation error $\left|\hat{\mathsf{w}}_{\{i,j\}}-w_{\{i,j\}}\right|$ is uniformly bounded by $\gamma/2$, an acyclic $G$ can be recovered from $\{\hat{\mathsf{w}}_{\{i,j\}}:\{i,j\}\in E\}$, similar to that from $\{w_{\{i,j\}}:\{i,j\}\in E\}$.

\begin{thm}\label{thm:Main}
Assume that $G=(V,I)$ is acyclic and for any $\{i,j\}\in I$ we have $|\beta_{\{i,j\}}|\in [\lambda,\mu]$ for some $\mu \geq \lambda >0$. Let $\mathsf{G}'=(V,\mathsf{T})$ be a maximum-weight spanning tree with respect to the weight assignment $\hat{\mathsf{w}}$. If
\begin{align}
\left|\hat{\mathsf{w}}_{\{i,j\}}-w_{\{i,j\}}\right|<\frac{\gamma}{2}, \quad \forall \{i,j\}\in E,\label{eq:assump}
\end{align}
then for any two distinct $i,j\in V$ we have $\{i,j\}\in I$ if and only if $\{i,j\}\in \mathsf{T}$ and $\hat{\mathsf{w}}_{\{i,j\}}>\gamma/2$.
\end{thm}

\begin{proof}
By Proposition~\ref{prop:PC1}, we have $w_{\{i,j\}} \geq \gamma >0= w_{\{k,l\}}$ for any $\{i,j\}\in I$ and $\{k,l\}\in D$. Under assumption \eqref{eq:assump}, this implies that $\hat{\mathsf{w}}_{\{i,j\}} > \gamma/2>\hat{\mathsf{w}}_{\{k,l\}}$ for any $\{i,j\}\in I$ and $\{k,l\}\in D$. By Proposition~\ref{prop:PC2}, we have $w_{\{i_1,i_{m+1}\}}\leq w_{\{i_s,i_{s+1}\}}-\gamma$ for any path $\left\{\{i_1,i_2\},\{i_2,i_3\},\ldots,\{i_m,i_{m+1}\}\right\}$ of length $m \geq 2$ in $G$ and any $s\in\{1,\ldots,m\}$. Under assumption \eqref{eq:assump}, this implies that $\hat{\mathsf{w}}_{\{i_1,i_{m+1}\}}<\hat{\mathsf{w}}_{\{i_s,i_{s+1}\}}$ for any path $\left\{\{i_1,i_2\},\{i_2,i_3\},\ldots,\{i_m,i_{m+1}\}\right\}$ of length $m \geq 2$ in $G$ and any $s\in\{1,\ldots,m\}$. The theorem thus follows directly from Proposition~\ref{prop:Cond} with $u=\hat{\mathsf{w}}$ and $\eta=\gamma/2$.
\end{proof}

We then establish the following algorithm to detect $G$ based on Theorem~\ref{thm:Main}.
\begin{table}[h]
\begin{center}
\begin{tabular}{ lc }
\hline\hline
\textbf{Algorithm 1} \\
\hline
\noindent \textbf{Input:} $(\mathsf{x}^d[t],\mathsf{y}[t])$, $t=1,\ldots,n$ and $(\mu,\lambda)$ such that $\mu \geq \lambda >0$.\\
\noindent \textbf{Output:} $\hat{\mathsf{G}}=(V,\hat{\mathsf{I}})$.\\
1 \ For all $1 \leq i < j \leq d$, compute $\hat{\mathsf{w}}_{\{i,j\}}$ according to \eqref{eq:What}.\\
\ \ \ \hspace{-2pt} Compute $\gamma$ from $(\mu,\lambda)$ according to \eqref{eq:gamma}.\\
2 \ Find a maximum-weight spanning tree $\hat{\mathsf{G}}'=(V,\hat{\mathsf{T}})$ over \\
\ \ \ \hspace{2pt}the vertex set $V$ with respect to the weight assignment \\
\ \ \ \hspace{2pt}$(\hat{\mathsf{w}}_{\{i,j\}}:1 \leq i < j \leq d)$.\\
3 \ Return $\hat{\mathsf{G}}=(V,\hat{\mathsf{I}})$ with $\hat{\mathsf{I}}=\left\{\{i,j\}\in \hat{\mathsf{T}}: \hat{\mathsf{w}}_{\{i,j\}}> \gamma/2\right\}.$\\
\hline\hline
\end{tabular}
\end{center}
\end{table}

The sample complexity of the algorithm is summarized in the following theorem.
\begin{thm}
Assume that $G=(V,I)$ is acyclic and for any $\{i,j\}\in I$ we have $|\beta_{i,j}|\in [\lambda,\mu]$ for some $\mu \geq \lambda >0$. Fix $0<\epsilon<1$ and let $n$ be a positive integer such that
\begin{align}
n & \geq \frac{128}{\gamma^2}\log\frac{d^2}{\epsilon}= \frac{64\pi d}{\left[\sigma(\lambda+3\mu)-\sigma(-\lambda+3\mu)\right]^2}\log\frac{d^2}{\epsilon}\label{samplecomplexity}.
\end{align}
Then with probability at least $1-\epsilon$, the algorithm can successfully detect the graph $G$ from $n$ i.i.d. samples of $(\mathsf{x}^d,\mathsf{y})$.
\end{thm}

\begin{proof}
By Proposition~\ref{prop:Chernoff} and Theorem~\ref{thm:Main}, we have
\begin{align*}
\mathrm{Pr}\left(\hat{\mathsf{G}} = G\right) & \geq \mathrm{Pr}\bigg(\bigcap_{\{i,j\}\in E}\Big\{\left|\hat{\mathsf{w}}_{i,j}-w_{i,j}\right|<\frac{\gamma}{2}\Big\}\bigg)\\
& \geq 1-d(d-1)e^{-\frac{n\gamma^2}{128}} \geq 1-d^2\cdot \frac{\epsilon}{d^2}=1-\epsilon.
\end{align*}
This completes the proof of the theorem.
\end{proof}

\section{Extension to Models with both Individual Effects and Cooperative Interactions}\label{sec:Ext}
In this section, we extend the results of Sections~\ref{sec:Id} and \ref{sec:Det} to the models that include both individual effects and cooperative interactions. More specifically, we shall assume that the conditional probability of the outcome $\mathsf{y}$ given the covariates $\mathsf{x}_V$ are given by:
\begin{align}
\mathrm{Pr}\left(\mathsf{y}=+1| \mathsf{x}_V=x_V\right) & =\sigma\Big(\sum_{i \in V}\beta_i x_i+\sum_{\{i,j\}\in E}\beta_{\{i,j\}}x_ix_j\Big),\\
\mathrm{Pr}\left(\mathsf{y}=-1| \mathsf{x}_V=x_V\right) & =\sigma\Big(-\sum_{i \in V}\beta_i x_i-\sum_{\{i,j\}\in E}\beta_{\{i,j\}}x_ix_j\Big)
\end{align}
for some real constants $\beta_i$'s and $\beta_{\{i,j\}}$'s. For any $i\in V$, we say that the covariate $\mathsf{x}_i$ has an {\em individual effect} on $\mathsf{y}$ if and only if $\beta_i\neq 0$; for any $\{i,j\}\in E$, we say that the covariates $\mathsf{x}_i$ and $\mathsf{x}_j$ {\em interact} if and only if $\beta_{\{i,j\}} \neq 0$. Let $\tilde{V}:=V\cup\{0\}$, $\tilde{E}:=\{\{i,j\}: i,j\in\tilde{V},i \neq j\}$, and $\beta_{\{i,0\}}:=\beta_i$ for all $i\in V$. Then, the structure of the model (including both individual effects and cooperative interactions) is fully captured by the simple graph $\tilde{G}=(\tilde{V},\tilde{I})$, where $\tilde{I}:=\{(i,j)\in \tilde{E}:\beta_{\{i,j\}}\neq 0\}$. As before, our goal is to recover the graph $\tilde{G}$ from i.i.d. samples of $(\mathsf{x}_V,\mathsf{y})$.

Toward this goal, we shall introduce an additional covariate $\mathsf{x}_0$, which we assume to be uniformly over $\{+1,-1\}$ and independent of $\mathsf{x}_V$, and use it to define an {\em auxiliary} model, for which the conditional probabilities of the outcome $\tilde{\mathsf{y}}$ given the covariates $\mathsf{x}_{\tilde{V}}$ are given by:
\begin{align}
\mathrm{Pr}\left(\tilde{\mathsf{y}}=+1| \mathsf{x}_{\tilde{V}}=x_{\tilde{V}}\right) & =\sigma\Big(\sum_{\{i,j\}\in \tilde{E}}\beta_{\{i,j\}}x_ix_j\Big),
\\
\mathrm{Pr}\left(\tilde{\mathsf{y}}=-1| \mathsf{x}_{\tilde{V}}=x_{\tilde{V}}\right) & =\sigma\Big(-\sum_{\{i,j\}\in \tilde{E}}\beta_{\{i,j\}}x_ix_j\Big).
\end{align}
By the results of Section~\ref{sec:Id}, if the underlying graph $\tilde{G}$ is known to be acyclic, it can be recovered from the weight assignment:
\begin{align}
\tilde{w}_{\{i,j\}}:=\left|2\mathrm{Pr}\left(\tilde{\mathsf{y}}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1\right)-1\right|
\end{align}
for all $\{i,j\}\in \tilde{E}$.

Note that when $j=0$, we trivially have
\begin{align*}
\mathrm{Pr}&\left(\tilde{\mathsf{y}}=+1|\mathsf{x}_i=+1,\mathsf{x}_0=+1\right)=\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_i=+1\right),
\end{align*}
so
\begin{align}
\tilde{w}_{\{i,0\}}=\left|2\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_i=+1\right)-1\right|\label{eq:WW1}
\end{align}
for all $i \in V$. On the other hand, when $i,j \neq 0$ and $i \neq j$, we may write
\begin{align*}
\mathrm{Pr}&\left(\tilde{\mathsf{y}}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1\right)\\
&=\frac{1}{2}\mathrm{Pr}\left(\mathsf{\tilde{y}}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1,\mathsf{x}_0=+1\right)+\\
&\hspace{40pt}\frac{1}{2}\mathrm{Pr}\left(\mathsf{\tilde{y}}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1,\mathsf{x}_0=-1\right)\\
&=\frac{1}{2}\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1\right)+\\
&\hspace{40pt}\frac{1}{2}\mathrm{Pr}\left(\tilde{\mathsf{y}}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1,\mathsf{x}_0=-1\right).
\end{align*}
To proceed, further note that
\begin{align}
&2^{d-2}\mathrm{Pr}\left(\tilde{\mathsf{y}}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1,\mathsf{x}_0=-1\right)\nonumber\\
&= \sum_{x_V:(x_i,x_j)=(+1,+1)}\mathrm{Pr}\left(\tilde{\mathsf{y}}=+1|\mathsf{x}_V=x_V,\mathsf{x}_0=-1\right)\nonumber\\
& = \sum_{x_V:(x_i,x_j)=(+1,+1)}\sigma\Big(-\sum_{k \in V}\beta_k x_k+\sum_{\{k,l\}\in E}\beta_{\{k,l\}}x_kx_l\Big)\nonumber\\
& = \sum_{x_V:(x_i,x_j)=(-1,-1)}\sigma\Big(\sum_{k \in V}\beta_k x_k+\sum_{\{k,l\}\in E}\beta_{\{k,l\}}x_kx_l\Big)\label{eq:nC}\\
&= \sum_{x_V:(x_i,x_j)=(-1,-1)}\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_V=x_V\right)\nonumber\\
&= 2^{d-2}\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_i=-1,\mathsf{x}_j=-1\right)\nonumber
\end{align}
where \eqref{eq:nC} follows from the simple change of variable $x_V \rightarrow -x_V$. It thus follows that
\begin{align*}
\mathrm{Pr}&\left(\tilde{\mathsf{y}}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1,\mathsf{x}_0=-1\right)\\
&\hspace{60pt}=\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_i=-1,\mathsf{x}_j=-1\right)
\end{align*}
and hence
\begin{align*}
\mathrm{Pr}&\left(\tilde{\mathsf{y}}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1\right)\\
&=\frac{1}{2}\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1\right)+\\
& \hspace{40pt}\frac{1}{2}\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_i=-1,\mathsf{x}_j=-1\right)
\end{align*}
giving
\begin{align}
\tilde{w}_{\{i,j\}}&=\left|\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_i=+1,\mathsf{x}_j=+1\right)+\right.\nonumber\\
&\hspace{40pt} \left.\mathrm{Pr}\left(\mathsf{y}=+1|\mathsf{x}_i=-1,\mathsf{x}_j=-1\right)-1\right|\label{eq:WW2}
\end{align}
for any $\{i,j\}\in E$.

Combining \eqref{eq:WW1} and \eqref{eq:WW2}, we conclude that the weight assignment $(\tilde{w}_{\{i,j\}}:\{i,j\}\in \tilde{E})$, and hence any acyclic $\tilde{G}$, can be fully recovered from the low-order joint probabilities $(p(\mathsf{x}_i,\mathsf{x}_j,\mathsf{y}):\{i,j\}\in E)$. The results of Section~\ref{sec:Det} can be extended similarly; the details are omitted.

\section{Simulation Results}\label{sec:Sim}
In our simulations, we randomly generate 5,000 logistic regression models, each including 10 independent binary covariates. For each model that we generate, we randomly choose 5 individual effects and 5 interaction pairs, resulting in an acyclic graph as its underlying structure. The model parameters $\beta_i$ for each individual effect and $\beta_{\{i,j\}}$ for each interaction pair are randomly assigned from a uniform distribution over $[-\mu,-\lambda]\cup[\lambda,\mu]$ with $0< \lambda<\mu$. In Fig.~\ref{fig:parameterranges}, we compare the detection rate of Algorithm~1 for different parameter ranges $(\lambda,\mu)=(0.3,0.5)$, $(0.5,1)$, and $(1,2)$, under the sample sizes of 300, 600, 900, 1,200 and 1,500, respectively. Here, we emphasize that a logistic regression model is correctly detected if and only if all 10 features, including 5 individual effects and 5 interaction pairs, are correctly detected. Clearly, the detection rate increases as the lower and upper bounds of the parameter range increase, and Algorithm~1 can achieve a high detection rate (at least 93\%) with a reasonably large number of training samples (around 1200). Also in~Fig.~\ref{fig:parameterranges}, we plot the fitted curves of the detection rate with respect to the increasing sample size based on the functional form $n \propto \log (1/(1-\textrm{detection rate}))$ that we derived in~(\ref{samplecomplexity}) in Section~\ref{sec:Det}. It is clear that the curves fit very well with the empirical results, thus validating the order-tightness of the lower bound in~(\ref{samplecomplexity}).

Next, we shall compare the performance of our algorithm (Algorithm~1) with three generic feature selection \mbox{algorithms:}  mRMR forward selection \cite{pe2005}, mutual information (MI) ranking \cite{sa2007}, and the problem-specific $L_1$-regularized logistic regression algorithm \cite{leesuin2006,parkthastie2007}. For the \mbox{mRMR} forward selection, MI ranking, and $L_1$-regularized logistic regression algorithms, we shall view each of the single variables $\mathsf{x}_i$'s and the interaction terms $\mathsf{x}_i\mathsf{x}_j$'s as a separate feature, and assume that the number of features to be selected is known. The estimation of mutual information in the mRMR forward selection and MI ranking algorithms is based on \cite{ant2001, pa2003}. For the $L_1$-regularized logistic regression algorithm, we tune the regularization parameter till the desired number of nonzero coefficients is obtained. 

\begin{figure}
\centering
\includegraphics[width=0.42\textwidth]{parameterrangeswithcurves}
\caption{The detection rates of Algorithm~1 for different parameter ranges $[-\mu,-\lambda]\cup[\lambda,\mu]$.}
\label{fig:parameterranges}
\end{figure}

In Fig.~\ref{fig:fourmethods-detection}, we compare the detection rate of Algorithm~1 with that of mRMR forward selection, MI ranking, and $L_1$-regularized logistic regression algorithms for $(\lambda,\mu)=(0.3,0.5)$ under a finite number of data samples. As we can see, Algorithm~1 achieves a significantly higher detection rate than the mRMR forward selection and the MI ranking algorithms, especially when the sample size is relatively small. This is due to the facts that: 1) Algorithm 1 exploits the fact that the underlying interaction graph is acyclic while the other two algorithms do not; 2) the proposed influence measure is much easier to estimate than MI. (The performances of the mRMR forward selection and the MI ranking algorithms are nearly identical since the feature candidates are pairwise independent.) On the other hand, the performances of Algorithm~1 and the $L_1$-regularized logistic regression algorithm appear to be very comparable. It is somewhat surprising that the $L_1$-regularized logistic regression algorithm performs well for {\em typical} problem instances with finite sample sizes. Intuitively, this is related to the fact that acyclic graphs are ``sparse" graphs. However, analyzing the sample complexity of the $L_1$-regularized logistic regression algorithm appears to be very challenging.

\begin{figure}
\centering
\includegraphics[width=0.42\textwidth]{fourmethods-detection}
\caption{Comparison of the detection rates among mRMR forward selection, MI ranking, $L_1$-regularized logistic regression, and Algorithm~1.}
\label{fig:fourmethods-detection}
\end{figure}


For completeness, in Fig.~\ref{fig:fourmethods-falsepositive} and Fig.~\ref{fig:fourmethods-prediction} we also compare the miss detection/false positive rate and prediction accuracy of \mbox{Algorithm~1} with those of the mRMR forward selection, the MI ranking, and the $L_1$-regularized logistic regression algorithms. Note that since the number of features selected is fixed in this case, the miss detection and false positive rates are identical. The prediction accuracy is calculated as follows. For each logistic regression model that we generate, we first obtain the model structure via one of the algorithms under the consideration, followed by standard logistic regression for parameter estimation. Once each logistic regression model is reconstructed, we randomly generate 200 testing samples to estimate the accuracy of the {\em outcome} prediction. As we can see, the relative performance among these algorithms is very similar to the case with the detection rate.




\begin{figure}
\centering
\includegraphics[width=0.42\textwidth]{fourmethods-falsepositive}
\caption{Comparison of the false positive rates for detecting nonzero model parameters among mRMR forward selection, MI ranking, $L_1$-regularized logistic regression, and Algorithm~1.}
\label{fig:fourmethods-falsepositive}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.42\textwidth]{fourmethods-prediction}
\caption{Comparison of the prediction accuracy among mRMR forward selection, MI ranking, $L_1$-regularized logistic regression, and Algorithm~1.}
\label{fig:fourmethods-prediction}
\end{figure}


\section{Concluding Remarks}\label{sec:Con}
An important problem in bioinformatics is to identify interactive effects among profiled variables for outcome prediction. In this paper, a simple logistic regression model with pairwise interactions among the binary covariates was considered. Modeling the structure of the interactions by a graph $G$, our goal was to recover $G$ from i.i.d. samples of the covariates and the outcome. When viewed as a feature selection problem, a simple quantity called influence is proposed as a measure of the marginal effects of the interaction terms on the outcome. For the case where $G$ is known to be acyclic, it is shown that a simple algorithm that is based on a maximum-weight spanning tree with respect to the plug-in estimates of the influences not only has strong theoretical performance guarantees, but can also outperform generic feature selection algorithms for recovering the graph from i.i.d. samples of the covariates and the outcome. A sample complexity analysis for detecting the interaction graph was provided, and the results were further extended to the model that includes both individual effects and pairwise interactions.

In our future work, we would like to understand the behavior of the $L_1$-regularized logistic regression algorithm from a theoretical standpoint, and also extend our results to the more challenging case where the interaction graph might be cyclic.

