Assessing Statistical Reliability of LiNGAM via Multiscale Bootstrap
Structural equation models and Bayesian networks have been widely used to study causal relationships between continuous variables. Recently, a non-Gaussian method called LiNGAM was proposed to discover such causal models and has been extended in various directions. An important problem with LiNGAM is that the results are affected by the random sampling of the data as with any statistical method. Thus, some analysis of the confidence levels should be conducted. A common method to evaluate a confidence level is a bootstrap method. However, a confidence level computed by ordinary bootstrap is known to be biased as a probability-value ($p$-value) of hypothesis testing. In this paper, we propose a new procedure to apply an advanced bootstrap method called multiscale bootstrap to compute $p$-values of LiNGAM outputs. The multiscale bootstrap method gives unbiased $p$-values with asymptotic much higher accuracy. Experiments on artificial data demonstrate the utility of our approach.
MBlingam_arXiv.tex
\documentclass{article}

\usepackage{makeidx}  % allows for indexgeneration
\usepackage{graphicx}
\usepackage{cite}

\usepackage{color}
\usepackage{amsmath}
\usepackage{epsfig}
\usepackage{url}

\usepackage{bm}

\newtheorem{thm}{Theorem}
\newtheorem{asp}{Assumption}
\newtheorem{lem}{Lemma}
\newtheorem{prp}{Proposition}
\newtheorem{proof}{Proof}
\renewcommand{\theproof}{}

%% Definitions of handy macros can go here

%\usepackage{graphics}
%\usepackage{color}
%\usepackage{amsmath}
%\usepackage{epsfig}
%\usepackage{url}

%%%%%%%%%%%%%%%%%%%%%%
%\usepackage[dvipdfm,bookmarks=true,bookmarksnumbered=true,%
%bookmarkstype=toc]{hyperref}
%\AtBeginDvi{\special{pdf:tounicode EUC-UCS2}}
%%%%%%%%%%%%%%%%%%%%%%

\newcommand{\qed}{\fbox{\rule[1pt]{0pt}{0pt}}}
\newcommand{\bsquare}{\hbox{\rule{6pt}{6pt}}}

\newcommand{\dataset}{{\cal D}}
\newcommand{\fracpartial}[2]{\frac{\partial #1}{\partial  #2}}

\newcommand{\bfa}{{\bf a}}
\newcommand{\bfd}{{\bf d}}
\newcommand{\bfe}{{\bf e}}
\newcommand{\bfx}{{\bf x}}
\newcommand{\bfy}{{\bf y}}
\newcommand{\bfz}{{\bf z}}
\newcommand{\bfr}{{\bf r}}

\def\bmu{\mbox{\boldmath $\mu$}}
\def\bsigma{\mbox{\boldmath $\sigma$}}
\def\bSigma{\mbox{\boldmath $\Sigma$}}
\def\bTheta{\mbox{\boldmath $\Theta$}}
\def\btheta{\mbox{\boldmath $\theta$}}
\def\bpsi{\mbox{\boldmath $\psi$}}
\def\bfeta{\mbox{\boldmath $\eta$}}
\def\bla{\mbox{\boldmath $\lambda$}}
\def\bgamma{\mbox{\boldmath $\gamma$}}
\def\bnu{\mbox{\boldmath $\nu$}}
\def\bal{\mbox{\boldmath $\alpha$}}
\def\bep{\mbox{\boldmath $\epsilon$}}
\def\bzeta{\mbox{\boldmath $\zeta$}}
\def\btau{\mbox{\boldmath $\tau$}}
\def\bbeta{\mbox{\boldmath $\beta$}}
\def\bphi{\mbox{\boldmath $\phi$}}
\def\bxi{\mbox{\boldmath $\xi$}}
\def\bupsilon{\mbox{\boldmath $\upsilon$}}
\def\bomega{\mbox{\boldmath $\omega$}}
\def\bpi{\mbox{\boldmath $\pi$}}
\def\bkappa{\mbox{\boldmath $\kappa$}}

\newcommand{\Wp}{{\widetilde{\bf W}}}
\newcommand{\perm}{{\bf P}}
\newcommand{\Be}{{\widehat{\bf B}}}
\newcommand{\Bp}{{\widetilde{\bf B}}}

\newcommand{\R}{{\bf R}}
\newcommand{\A}{{\bf A}}
\newcommand{\B}{{\bf B}}
\newcommand{\C}{{\bf C}}
\newcommand{\E}{{\bf E}}
\newcommand{\M}{{\bf M}}
\newcommand{\J}{{\bf J}}
\newcommand{\K}{{\bf K}}
\newcommand{\W}{{\bf W}}
\newcommand{\I}{{\bf I}}
\newcommand{\D}{{\bf D}}
\newcommand{\G}{{\bf G}}
\newcommand{\x}{{\bf x}}
\newcommand{\X}{{\bf X}}
\newcommand{\bP}{{\bf P}}
\newcommand{\bS}{{\bf S}}
\newcommand{\bfL}{{\bf L}}
\newcommand{\T}{{\bf T}}
\newcommand{\U}{{\bf U}}
\newcommand{\V}{{\bf V}}
\newcommand{\y}{{\bf y}}
\newcommand{\Y}{{\bf Y}}
\newcommand{\Z}{{\bf Z}}
\newcommand{\Q}{{\bf Q}}
\newcommand{\s}{{\bf s}}
\newcommand{\F}{{\bf F}}
\newcommand{\bLambda}{{\bf \Lambda}}

%\title{Assessing Statistical Reliability of LiNGAM via Multiscale Bootstrap}
\title{Computing $p$-values of LiNGAM outputs via Multiscale Bootstrap}
\author{Yusuke Komatsu\thanks{Department of Mathematical and Computing Sciences, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8552, Japan.}, Shohei Shimizu\thanks{The Institute of Scientific and Industrial Research (ISIR), Osaka University, Mihogaoka 8-1, Ibaraki, Osaka 567-0047, Japan. Email: sshimizu@ar.sanken.osaka-u.ac.jp}, Hidetoshi Shimodaira\thanks{Department of Mathematical and Computing Sciences, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8552, Japan}
}
\date{}


% The \author macro works with any number of authors. There are two commands
% used to separate the names and addresses of multiple authors: \And and \AND.
%
% Using \And between authors leaves it to \LaTeX{} to determine where to break
% the lines. Using \AND forces a linebreak at that point. So, if \LaTeX{}
% puts 3 of 4 authors names on the first line, and the last on the second
% line, try using \AND instead of \And before the third author name.

\newcommand{\fix}{\marginpar{FIX}}
\newcommand{\new}{\marginpar{NEW}}

\begin{document}

\maketitle

\begin{abstract}
Structural equation models and Bayesian networks have been widely used to study causal relationships between continuous variables.
Recently, a non-Gaussian method called LiNGAM was proposed to discover such causal models and has been extended in various directions. An important problem with LiNGAM is that the results are affected by the random sampling of the data as with any statistical method. Thus, some analysis of the confidence levels should be conducted. A common method to evaluate a confidence level is a bootstrap method. However, a confidence level computed by ordinary bootstrap is known to be biased as a probability-value ($p$-value) of hypothesis testing. In this paper, we propose a new procedure to apply an advanced bootstrap method called multiscale bootstrap to compute $p$-values of LiNGAM outputs. The multiscale bootstrap method gives unbiased $p$-values with asymptotic much higher accuracy. Experiments on artificial data demonstrate the utility of our approach.
 \end{abstract}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%An important goal in many empirical sciences is to understand causal mechanisms underlying their objective systems such as human social behavior and natural phenomena. 
Structural equation models \cite{Bollen89book} and Bayesian networks \cite{Pearl00book,Spirtes93book} have been widely applied to analyze causal relationships in many fields. 
Many methods \cite{Pearl00book,Spirtes93book} have been developed to discover such a causal model when no prior knowledge on the network structure is available. 
Recently, a non-Gaussian method called LiNGAM \cite{Shimizu06JMLR} was proposed. 
The new method estimates a causal ordering of variables using passive observational data alone. 
The estimated ordering is correct if the causal relations form a linear structural equation model with non-Gaussian external influence variables {\it and} the sample size is {\it infinitely large}. 
In practice, however, the sample size is finite. 
The finite sample size induces statistical errors in the estimation, and the estimated ordering may not be right even when the model assumptions are reasonable. 
Thus, some analysis of the statistical reliability or confidence level of the estimated ordering should be done. 
%In other words, we want to know if the estimated ordering and its suborderings are reliable, {\it i.e.}, likely to be found in most of datasets created by random sampling from the same model.  
In this paper, we discuss such reliability analysis of LiNGAM. 
%Most of more recent works \cite{Hoyer09NIPS,Lacerda08UAI} on LiNGAM have focused on identifiability to generalize the model and paid little attention to testing and reliability, although all of the issues are important to apply LiNGAM on real-world data. 

A common procedure to evaluate such a confidence level is statistical hypothesis testing \cite{Lehmann08Test}. 
In statistical testing, one computes a probability-value ($p$-value) of a hypothesis. The hypothesis is rejected when the $p$-value is  not greater than a pre-specified level of significance, say 5\%. 
%The hypothesis is more reliable for the larger the $p$-value. 
There are several approaches to define a $p$-value. 
Bootstrapping \cite{Efron93book} is a well-known computational method for computing confidence levels when a simple mathematical formula is difficult to derive. 
It is a resampling method to approximate a random sample by a bootstrap sample that is created by random sampling with replacement from the original single dataset. 
Felsenstein \cite{Felsenstein85Evol} proposed to use bootstrapping to define a $p$-value in the context of phylogenetic tree selection of molecular evolution in bioinformatics. 
He defined a $p$-value of a tree by a frequency called bootstrap probability that the tree is found to be optimal when tree selection is performed for a number of bootstrap replicates of the original dataset. 
The idea has been applied to other multivariate analyses including Bayesian networks  \cite{Fried99UAI}. % and cluster analysis \cite{Suzuki06PV}. 

However, it is known that the bootstrap probability is biased as a $p$-value \cite{Efron96boot,Hillis93SB}. 
The naive bootstrapping tends to give overconfidence in wrong hypotheses. % or trees.  
Thus, some advanced bootstrap methods to achieve higher accuracy have been proposed \cite{Efron96boot,Hall92book,Shimo02SB,Shimo08JSPI}. 
Among others, multiscale bootstrapping \cite{Shimo02SB,Shimo08JSPI} is much more accurate but still easy to implement and has been successful in the field of phylogenetic tree selection. 

In this paper, we propose to apply the multiscale bootstrap to compute confidence levels, {\it i.e.}, $p$-values, of variable orderings estimated by LiNGAM. 
The paper is structured as follows. 
First, in Section~\ref{sec:background}, we briefly review LiNGAM and multiscale bootstrap. 
In Section~\ref{sec:MB-lingam} we propose a new procedure to compute $p$-values of LiNGAM outputs using the multiscale bootstrap method. % and in Section~\ref{sec:relatedwork} discuss related work.  
The multiscale bootstrap method is tested using artificial data in Section~\ref{sec:simulations}. 
Conclusions are given in Section~\ref{sec:conc}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Background}\label{sec:background}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{LiNGAM}\label{sec:lingam}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In \cite{Shimizu06JMLR}, a non-Gaussian variant of structural equation models and Bayesian networks, which is called LiNGAM, was proposed. 
Assume that observed data are generated from a process represented graphically by a directed acyclic graph, {\it i.e.}, DAG. 
Let us represent this DAG by a $m$$\times$$m$  adjacency matrix $\B$$=$$\{b_{ij}\}$ where every $b_{ij}$ represents the connection strength from a variable $x_j$ to another $x_i$ in the DAG, {\it i.e.}, the {\it direct} causal effect of $x_j$ on $x_i$. 
Let us further define $\A=(\I-\B)^{-1}$. The $(j,i)$-element $a_{ji}$ represents the {\it total} causal effect of $x_i$ on $x_j$ \cite{Hoyer07IJAR}.    
Moreover, let us denote by $k(i)$ a causal order of variables $x_i$ in the DAG so that no later variable influences any earlier variable. 
For example, a variable $x_j$ is not causally influenced by a variable $x_i$, {\it i.e.}, $a_{ji}$$=$$0$, if $k(j)<k(i)$.
Moreover, assume that the relations between variables are linear. 
%Without loss of generality, each observed variable $x_i$ is assumed to have zero mean. 
Then we have
\begin{equation}
x_i = \sum_{k(j)<k(i)} b_{ij}x_j + e_i,\label{eq:model0}
\end{equation}
where $e_i$ is an external influence variable. 
All external influences $e_i$ are continuous random variables having \emph{non-Gaussian} distributions with zero means and non-zero variances, and $e_i$ are independent of each other so that there is no unobserved confounding variables \cite{Spirtes93book}. 
We emphasize that $k(j)$$<$$k(i)$ does not necessarily imply that $x_j$ influences $x_i$. 
It only implies that $a_{ji}$$=$$0$, and $a_{ij}$ can be either zero or non-zero. 
The causal ordering $k(i)$ only defines a {\it partial} order of variables, which is enough to define a DAG. 
In \cite{Shimizu06JMLR}, the LiNGAM model (\ref{eq:model0}) was shown to be identifiable without using any prior knowledge on the network structure. 
That is, the variable orders $k(i)$ and connection strengths $b_{ij}$ are estimable solely based on the data matrix of $\x=[x_1,\cdots,x_m]^T$ . 
%We rewrite the model (\ref{eq:model0}) in a matrix form as follows: 
%\begin{equation}
%\bfx = \B\bfx + \bfe,\label{eq:model}
%\end{equation} 
%where $\bfx$ is a $m$-dimensional variable vector, and $\B$ can be permuted by simultaneous equal row and column permutations to be strictly lower triangular due to the acyclicity assumption \cite{Bollen89book}. Strict lower triangularity is here defined as to have a lower triangular structure with all zeros on the diagonal. 
In \cite{Shimizu06JMLR}, a discovery algorithm based on independent component analysis (ICA) \cite{Hyva01book}, which is called LiNGAM algorithm, was also proposed to estimate $k(i)$ and $b_{ij}$. 
%In what follows, we refer to the ICA-based algorithm as LiNGAM algorithm. 
%Give a short explanation of the original LiNGAM algorithm. Do ICA + matrix permutation.

\subsection{Bootstrap probability}
Denote by $\x$ a $m$-dimensional random variable vector and by $\X$$=$$(\x_1,\cdots,\x_n)$ a random sample of size $n$ from the distribution of $\x$. 
Further, define a function $f(\X)$ so that $f(\X)$$=$$0$ if a hypothesis is rejected and otherwise $f(\X)$$=$$1$. 
Suppose that we obtain a $m$$\times$$n$ data matrix $\overline{\X}$ that is generated from $\x$, and the function $f(\overline{\X})$$=$$1$. 
Then, it is useful to evaluate how statistically reliable the value of $f(\overline{\X})$$=$$1$ is since the function could return 0 for another data matrix due to sample fluctuation. 
In \cite{Felsenstein85Evol}, Felesenstein proposed to use bootstrapping \cite{Efron93book} to evaluate such reliability. 
Let us denote by $\X_q^{\ast}$ a $q$-th bootstrap sample of size $n^{\ast}$, which is created by random sampling with replacement from the columns of $\X$. 
In ordinary bootstrap, $n^{\ast}$ is taken to be $n$. %, whereas in multiscale bootstrapping, a number of different values of $n^{\ast}$ are used.
Then, the bootstrap probability $p^{BP}$ is defined as a frequency that $f(\X^{\ast})$$=$1: 
\begin{eqnarray}
p^{BP} &=& \frac{1}{Q} \sum_{q=1}^{Q} f\left(\X_q^{\ast}\right),\label{eq:bp}
\end{eqnarray}
where $Q$ is the number of bootstrap replications. 
A testing procedure was proposed that the hypothesis is rejected if $p^{BP}$ is not greater than a significance level $\alpha$ ($0$$<$$\alpha$$<$$1$), say 0.05. 
However, it is known that $p^{BP}$ is {\it biased} as a $p$-value \cite{Efron96boot,Hillis93SB}. The multiscale bootstrap \cite{Shimo02SB,Shimo08JSPI} corrects the bias and gives a more accurate $p$-value. This is explained in more detail in the next subsection. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Unbiasedness}\label{sec:ub}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
To discuss the bias of a $p$-value, it is conventional \cite{Efron96boot} to assume that 
there {\it exists} a function $g$ that transforms a random sample $\X$ to a $K$-dimensional random vector $\y$$=$$[y_1,\cdots,y_K]^T$ that (at least approximately) follows a Gaussian distribution with an unknown mean vector ${\bm \mu}$ and covariance identity $\I$, {\it i.e.}, $N_K({\bm \mu},\I)$.  
Note that it is {\it not} necessary to specify the actual functional form of $g$ and dimension $K$. 
Let us denote by $\mathcal{H}$ such a class of $\y$ that $f(\X)$$=$$1$.  
Then, the null hypothesis $f(\X)$$=$$1$ can be described as ${\bm \mu}$$\in$$\mathcal{H}$ in terms of a region in the parameter space. 
We only have to consider $\y$ to discuss the bias of a $p$-value computed based on $\X$ due to the {\it  transformation-respecting property} of bootstrapping \cite{Efron93book}. 


In statistical hypothesis testing, the null hypothesis ${\bm \mu}$$\in$${\mathcal H}$ is rejected when a $p$-value computed based on $\y$, which is denoted by $p(\y)$, is not greater than a significance level $\alpha$. 
A test controls a type-I error if the probability of false rejection under the null hypothesis is not greater than $\alpha$. 
This is a desirable property of a testing procedure. 
Another desirable property is {\it unbiasedness} \cite{Lehmann08Test}.  
A test is unbiased if the  probability of correct rejection under alternative hypotheses is not less than $\alpha$, and the type-I error is also controlled. 
Then an unbiased test is formally defined to be a test that uses a $p$-value $p(\y)$ satisfying
\begin{eqnarray}
{\rm Prob}\{p(\y)<\alpha\}\le \alpha, \hspace{1mm} {\bm \mu} \in \mathcal{H}\hspace{3mm}{\rm and}\hspace{3mm}{\rm Prob}\{p(\y)<\alpha\}\ge \alpha, \hspace{1mm} {\bm \mu} \notin \mathcal{H}.
\end{eqnarray}
%\begin{eqnarray}
%{\rm Prob}\{p(\y)<\alpha\}\le \alpha, \hspace{1mm} {\bm \mu} \in \mathcal{H}. \label{eq:3}\\
%{\rm Prob}\{p(\y)<\alpha\}\ge \alpha, \hspace{1mm} {\bm \mu} \notin \mathcal{H}.\label{eq:4}
%\end{eqnarray}
Let us denote by $\partial \mathcal{H}$ the boundary of $\mathcal{H}$. 
%To satisfy the inequalities (\ref{eq:3}) and (\ref{eq:4}), the following equation needs to hold \cite{Lehmann08Test}:
To satisfy the inequalities above, the following equation needs to hold \cite{Lehmann08Test}:
\begin{eqnarray}
{\rm Prob}\{p(\y)<\alpha\} = \alpha, \hspace{1mm} {\bm \mu} \in \partial \mathcal{H}.\label{eq:5}
\end{eqnarray}
In other words, $p(\y)$ follows a {\it uniform} distribution over the interval $[0,1]$. 
It has been shown \cite{Shimo02SB} that $p^{BP}$ has a rather large bias to meet the unbiasedness condition (\ref{eq:5}):
\begin{eqnarray}
{\rm Prob}\{p^{BP}(\y)<\alpha\} = \alpha + O(n^{-1/2}), 
\end{eqnarray}
where $O(\cdot)$ is the Landau symbol. 
%Multiscale bootstrap \cite{Shimo02SB} was developed to reduce the bias. 
Multiscale bootstrap \cite{Shimo02SB} reduces the bias. 
Let $p^{MB}$ denote a $p$-value computed by multiscale bootstrap. 
It can be shown that  $p^{MB}$ is approximately unbiased with asymptotic third-order accuracy:
\begin{eqnarray}
{\rm Prob}\{p^{MB}(\y)<\alpha\} = \alpha + O(n^{-3/2}),\label{eq:mb2}
\end{eqnarray}
Thus, multiscale bootstrap gives a $p$-value with much higher-order accuracy than ordinary bootstrap. %\footnote{The theory of multiscale bootstrap has been extended for exponential families \cite{Shimo04AS}.} 
Rigorously speaking, the boundary $\partial \mathcal{H}$ needs to be assumed to be smooth or approximately smooth. 
Otherwise, no unbiased test can be defined \cite{Lehmann08Test}. 
However, it has been shown that $p^{MB}$ is less biased than $p^{BP}$ even if the boundary is non-smooth \cite{Shimo08JSPI}. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Multiscale Bootstrap}\label{sec:MB}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In \cite{Shimo08JSPI}, the theory of multiscale bootstrap \cite{Shimo02SB} was extended, and a class of unbiased $p$-values including $p^{MB}$ in (\ref{eq:mb2}) was obtained. 
%Let us define a scale $\sigma$$=$$\sqrt{n/n^{\ast}}$. 
Let $\y^{\ast}$ denote the $\y$ vector corresponding to $\X^{\ast}$$=$$[\x_1^{\ast}, \cdots, \x_{n^{\ast}}^{\ast}]$. Then
the standard deviation of $\y^{\ast}$ is proportional to $1/\sqrt{n^{\ast}}$, and its value
relative to the case $n^{\ast}$$=$$n$ is called `scale' of bootstrap resampling; this is
defined by $\sigma$$=$$\sqrt{n/n^{\ast}}$.
Then the bootstrap probability $p^{BP}$ in (\ref{eq:bp}) is a function of $\sigma^2$, which is denoted by  $p^{BP}_{\sigma^2}$ for clarity. 
The fundamental idea of the extended multiscale bootstrap \cite{Shimo08JSPI} is to compute the bootstrap probability $p^{BP}_{\sigma^2}$ with the scale $\sigma^2$$=$$-1$, {\it i.e.}, the bootstrap sample size $n^{\ast}$$=$$-n$. 
Of course, it is impossible to set $n^{\ast}$$=$$-n$. 
Therefore, one first select several scales $\sigma$$>$$0$, computes the bootstrap probability for each of the corresponding bootstrap sample sizes $n/\sigma^2$ and extrapolates the bootstrap probabilities to  $\sigma^2$$=$$-1$, {\it i.e.}, $n^{\ast}$$=$$-n$. 

We now review a procedure to compute such unbiased $p$-values. 
Let us define a bootstrap $z$-value by 
\begin{eqnarray}
z_{\sigma^2}=-\Phi^{-1}\left(p^{BP}_{\sigma^2}\right),\label{eq:nz}
\end{eqnarray}
where $\Phi^{-1}$ is the inverse of the distribution function $\Phi$ of the standard Gaussian distribution $N(0,1)$. 
Further, let us call $\sigma z_{\sigma^2}$ a normalized bootstrap $z$-value. 
Then, consider to model the changes in $\sigma z_{\sigma^2}$ along the changing the scale $\sigma$ by a model $\psi(\sigma^2|{\bm \beta})$, where ${\bm \beta}$$=$$[\beta_0, \cdots, \beta_{h-1}]^T$ is a parameter vector of the model. 
Two model classes are proposed in \cite{Shimo08JSPI}:
\begin{eqnarray}
\psi_1^{h}(\sigma^2|{\bm \beta})  &=& \sum_{j=0}^{h-1} \beta_j \sigma^{2j}, \hspace{1mm} h\ge1. \label{eq:ploy}\\
\psi_2^{h}(\sigma^2|{\bm \beta})  &=& \beta_0 + \sum_{j=1}^{h-2} \frac{\beta_j\sigma^{2j}}{1+\beta_{h-1}(\sigma-1)}, \hspace{1mm}  0\le \beta_{h-1}\le1,\hspace{1mm}  h\ge3.\label{eq:sing}
\end{eqnarray}
The model (\ref{eq:ploy}) is reasonable when the boundary $\partial \mathcal{H}$ is smooth, and the model (\ref{eq:sing}) is preferable when $\mathcal{H}$ is a cone and $\partial \mathcal{H}$ is not smooth. 
To estimate the models, a number of sets of bootstrap replicates with different scales $\sigma_d$ ($d$$=$$1$, $\cdots$, $D$) are first created, and subsequently the bootstrap probability $p^{BP}_{\sigma_d^2}$ for each scale is computed. 
Note that the bootstrap sample sizes may be different from that of the original dataset. 
Then, a set of scales and normalized bootstrap $z$-values \{$\sigma_d$, $\sigma_d z_{\sigma_d^2}$\} is obtained. 
Note that $z_{\sigma_d^2}$ is computed based on $p^{BP}_{\sigma_d^2}$ using (\ref{eq:nz}). 
Finally, the model parameter vector ${\bm \beta}$ are estimated using the set of scales and normalized bootstrap $z$-values. 
The maximum likelihood method is applied since $Qp^{BP}_{\sigma^2}$ follows a binomial distribution. A best model  $\psi_{best}(\sigma^2|{\bm \beta})$ is selected using an information criterion AIC \cite{Akaike74AIC}.

Then, a class of $p$-values using the best model $\psi_{best}(\sigma^2|{\bm \beta})$ is derived:
\begin{eqnarray}
p^{MB}_h &=& \Phi\left\{-\sum_{j=0}^{h-1}\frac{(-1-\sigma^2_0)^j}{j!}\frac{\partial^j \psi_{best}(\sigma^2|\hat{{\bm \beta}})}{\partial (\sigma^2)^j}\left|_{\sigma_0^2}\right.\right\},\label{eq:mbk}
\end{eqnarray}
where $\sigma_0^2$ is taken to be unity. 
The right side of (\ref{eq:mbk}) is the first $h$ terms of the Taylor series of the slope of $z_{\sigma^2}$ at $1/\sigma=1$, {\it i.e.}, $\partial z_\sigma^2/\partial(1/\sigma)|_1$, around $\sigma_0^2$. 
It can be shown that $p^{MB}_2$ is actually equal to $p^{MB}$ in (\ref{eq:mb2}) that achieves the unbiasedness with asymptotic third-order accuracy. 
Further, $p^{MB}_1$ turns out to be the naive bootstrap probability $p^{BP}$ in (\ref{eq:bp}). 
The larger $h$ gives an unbiased $p$-value with asymptotic higher-order accuracy \cite{Shimo08JSPI}.
However, it also makes the maximum likelihood estimation less stable. 
In practice, $h$$=$$2$ or $3$ is often used. 

\begin{figure}[h]
\begin{center}
%\includegraphics[width=.5\textwidth,trim=0 350 0 0]{territory3.eps}
\includegraphics[width=2in]{region.eps}
\end{center}
\caption{Four regions $H_{12}^{+}$, $H_{12}^{-}$, $H_{21}^{+}$, and $H_{21}^{-}$.}
\label{fig:region}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{A multiscale bootstrap procedure to assessing reliability of LiNGAM}\label{sec:MB-lingam}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{Null hypotheses}\label{sec:MB-null}
%We first define null hypotheses to be tested in LiNGAM. 
We first define null hypotheses tested. 
We here focus on the following four types of hypotheses between $x_i$ and $x_j$ ($i$$\neq$$j$), although we can test hypotheses that describe the relations between more than two variables similarly: 
\begin{enumerate}
\item $H_{ij}^{+}$: a hypothesis that $x_i$ is directly caused by $x_j$, and its connection strength is positive, {\it i.e.}, $b_{ij}$$>$$0$; 
\item $H_{ij}^{-}$: a hypothesis that $x_i$ is directly caused by $x_j$, and its connection strength is negative, {\it i.e.}, $b_{ij}$$<$$0$; 
\item $H_{ji}^{+}$: a hypothesis that $x_j$ is directly caused by $x_i$, and its connection strength is positive, {\it i.e.}, $b_{ji}$$>$$0$; 
\item $H_{ji}^{-}$: a hypothesis that $x_j$ is directly caused by $x_i$, and its connection strength is negative, {\it i.e.}, $b_{ji}$$<$$0$. 
\end{enumerate}
%The parameter space is also divided to the four regions by the sign of the connection strength between $x_i$ and $x_j$ as well as their causal ordering. 
See Fig.~\ref{fig:region} for the four regions of the parameter space in {\it two} variable cases {\it around the origin} that the connection strengths $b_{12}$ and $b_{21}$ are zeros. 
LiNGAM outputs $k(2)$$>$$k(1)$ if $|b_{12}|$$<$ $|b_{21}|$, and otherwise $k(1)$$>$$k(2)$ since each total effect $a_{ij}$ is equal to the corresponding direct effect $b_{ij}$ in two variable cases \cite{Shimizu06JMLR}.  
We note that the signs of connections strengths are important and interesting in many applications  \cite{Bollen89book,Silva06JMLR} as well as the variable orderings.
Further, this way of dividing the space based on the signs and orderings would make the boundaries of the regions be closer to be smooth than solely based on the orderings and help the multiscale bootstrap work better. 
%In Fig.~\ref{fig:region}, each of the four regions $H_{12}^{+}$, $H_{12}^{-}$, $H_{21}^{+}$, and $H_{21}^{-}$ is a cone whose vertex is on the origin point and boundary is smooth other than the origin. 
%In contrast, each of unions  $H_{12}^{+}$$\cup$$H_{12}^{-}$ and $H_{21}^{+}$$\cup$$H_{21}^{-}$ is represented by a region whose boundary is much more complicated since it connects the corresponding two cone-shaped regions at the vertex (the origin). 
%(The issue of non-smoothness is not very specific to LiNGAM. A successful application of multiscale boostrap is tree selection \cite{Shimo02SB}, but hypotheses in tree selection also form non-smooth boundaries \cite{Shimo08JSPI}.)
%We experimentally confirmed the conjecture conducting small simulations. Future work would conduct a theoretical analysis of the conjecture. 

%\subsection{MB-LiNGAM procedure}
We now propose a new procedure to apply Multiscale Bootstrap to LiNGAM, which we call {\it MB-LiNGAM}:

\noindent
  %\pagebreak
  \rule{\columnwidth}{0.5mm}
       { \sffamily
	MB-LiNGAM procedure
	 \begin{enumerate}
	 \item \label{step:ica} Select the scales $\sigma_1$, $\cdots$, $\sigma_D$ ($D$$\ge$$2$) so that $n^{\ast}_d$$=$$n/\sigma_d^2$ is an integer and choose the number of bootstrap replicates $Q$. 
	 \item Generate $Q$ bootstrap replicates $\X^{\ast}_{q,d}$ ($q$$=$$1$, $\cdots$, $Q$) for each scale $\sigma_d$, {\it i.e.}, each bootstrap sample size $n^{\ast}_d$$=$$n/\sigma_d^2$.
	 \item\label{step:lingam} Perform LiNGAM algorithm to each bootstrap replicate $\X^{\ast}_{q,d}$ and compute the bootstrap probabilities $p^{BP}_d(H_{ij}^{+})$ and $p^{BP}_d(H_{ij}^{-})$ ($i$$\neq$$j$) for each scale $\sigma_d$, where $p^{BP}_d(H)$ denotes the bootstrap probability of a hypothesis $H$ for scale $\sigma_d$. 
	 \item Compute the multiscale bootstrap $p$-values $p^{MB}_h(H_{ij}^{+})$ and $p^{MB}_h(H_{ij}^{-})$ ($i$$\neq$$j$) using the procedure in Section \ref{sec:MB}, more specifically (\ref{eq:mbk}),  where $p^{MB}_h(H)$ denotes the multiscale bootstrap $p$-value of $H$ with the order $h$. %Here, $h$$=$$3$ is used. 
	 \end{enumerate}
       } \vspace{-4mm}
\noindent \rule{\columnwidth}{0.5mm}
In the simulations below, the ICA part of LiNGAM algorithm is run several times in Step~\ref{step:lingam}. 
Each time the initial point of the optimization is randomly changed. 
The set of the estimates that achieves the largest value of an ICA objective function is used in the subsequent steps. 
It is a common practice to alleviate the effects of possible local maxima. 


%\begin{figure}[htbp]
%\begin{tabular}{cc}
%\begin{minipage}{0.5\hsize}
%\begin{center}
%\includegraphics[width=2in]{region.eps}
%\caption{Four regions $H_{12}^{+}$, $H_{12}^{-}$, $H_{21}^{+}$, and $H_{21}^{-}$.}
%\label{fig:region}
%\end{center}
%\end{minipage}
%\begin{minipage}{0.5\hsize}
%\begin{center}
%\includegraphics[width=2in]{only6ss2.eps}
%\caption{Scatterplots of empirical rejection probabilities of $H_{32}^{+}$ by ordinary bootstrap and those by multiscale bootstrap versus significance levels in the six variable cases. 
%}
%\label{fig:plots}
%\end{center}
%\end{minipage}
%\end{tabular}
%\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Related work}\label{sec:relatedwork}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Related work}
Some methods have been proposed to test significance of direct effects $b_{ij}$\cite{Shimizu06JMLR,Hoyer07IJAR}.  
For simplicity, let us consider two variable cases, where each direct effect $b_{ij}$ is equal to the corresponding total effect $a_{ij}$ as mentioned above.
Those methods test if each of effects $b_{ij}$ is zero or not and imply that $k(i)$$<$$k(j)$ if `$b_{ij}$$=$$0$' is accepted and that $k(j)$$<$$k(i)$ if `$b_{ji}$$=$$0$' is accepted. Such a procedure would work if $b_{ij}$ or $b_{ji}$ is exactly zero. 
However, in reality, the assumptions of the model (\ref{eq:model0}) are more or less violated, and hence both of $b_{ij}$ and $b_{ji}$ could be non-zero. In such cases, those existing methods might reject both of the hypotheses and not give much information on which ordering is better. Even in the cases, our approach tells which ordering is better or statistically more reliable comparing bootstrap probabilities of the orderings. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Simulations}\label{sec:simulations}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We first created three LiNGAM models with $m$$=$$2$ variables: 
\begin{eqnarray}
\left[
\begin{array}{c}
x_1\\
x_2
\end{array}
\right]
& =& 
\left[
\begin{array}{cc}
0 & b\\
b & 0
\end{array}
\right]
\left[
\begin{array}{c}
x_1\\
x_2
\end{array}
\right]
+
\left[
\begin{array}{c}
e_1\\
e_2
\end{array}
\right],
\end{eqnarray}
where $b$$=$$0$, $0.01$ or $0.1$, and $e_1$ and $e_2$ followed a Laplace distribution with mean zero and variance two. 
%All the models are on boundaries of the hypothesis regions in the parameter space. 
The model with $b$$=$$0$ is on the boundary of $H_{12}^{+}$, $H_{12}^{-}$, $H_{21}^{+}$, and $H_{21}^{-}$. The model with $b$$=$$0.01$ or $0.1$ is on the boundary between $H_{12}^{+}$ and $H_{21}^{+}$. Further, we created two LiNGAM models with $m$$=$$6$ variables: 
\begin{eqnarray}
\left[
\begin{array}{c}
x_1\\
x_2\\
x_3\\
x_4\\
x_5\\
x_6
\end{array}
\right]
& =& 
\left[
\begin{array}{cccccc}
0 & 0 & 0 & 0 & 0 & 0\\
b & 0 & 0 & 0 & 0 & 0\\
b & 0 & 0 & 0 & 0 & 0\\
b & b & 0 & 0 & 0 & 0\\
0 & b & 0 & b & 0 & 0\\
b & b & b & 0 & b & 0
\end{array}
\right]
\left[
\begin{array}{c}
x_1\\
x_2\\
x_3\\
x_4\\
x_5\\
x_6
\end{array}
\right]
+
\left[
\begin{array}{c}
e_1\\
e_2\\
e_3\\
e_4\\
e_5\\
e_6
\end{array}
\right],
\end{eqnarray}
where $b$$=$$0$ or $0.5$, and $e_1$ and $e_2$ also followed a Laplace distribution whose mean zero and variance two. 
%All the models are on boundaries of the hypothesis regions as well. 
We randomly generated 1280 datasets with sample size 1000 under each of the five models. 
%each combination of number of variables and connection strength $b$. 
Then we applied MB-LiNGAM procedure in Section~\ref{sec:MB-lingam} to the datasets. 
The scales $\sigma_d$ were selected so that they gave integer values of bootstrap sample size and were (approximately) equally-spaced in log-scale between $1/9$ and $9$ ($d$$=$$1$, $\cdots$, $13$). 
The number of bootstrap replicates $Q$ was 1000, and the value $h$ for $p_h^{MB}$ was 3.

\begin{figure}[!b]
\centering
\begin{minipage}{.25\textwidth}
\centering
%\caption{
BP, $b=0$%}

\includegraphics[width=1.3in]{2epsfile/komatsu20090116a/rap1ss.eps}
\end{minipage}
\begin{minipage}{.25\textwidth}
\centering
%\caption{
BP, $b=0.01$%}

\includegraphics[width=1.3in]{2epsfile/komatsu20090116b/rap1bss.eps}
\end{minipage}
\begin{minipage}{.25\textwidth}
\centering
%\caption{
BP, $b=0.1$%}

\includegraphics[width=1.3in]{2epsfile/komatsu20090116c/rap1css.eps}
\end{minipage}

\vspace{5mm}

\begin{minipage}{.25\textwidth}
\centering
%\caption{
MB, $b=0$%}

\includegraphics[width=1.3in]{2epsfile/komatsu20090116a/mup1ss.eps}
\end{minipage}
\begin{minipage}{.25\textwidth}
\centering
%\caption{
MB, $b=0.01$%}

\includegraphics[width=1.3in]{2epsfile/komatsu20090116b/mup1bss.eps}
\end{minipage}
\begin{minipage}{.25\textwidth}
\centering
%\caption{
MB, $b=0.1$%}

\includegraphics[width=1.3in]{2epsfile/komatsu20090116c/mup1css.eps}
\end{minipage}
\caption{Top row: Histograms of $p$-values of $H_{21}^{+}$ by ordinary bootstrap (BP). Bottom row:  Histograms of $p$-values of $H_{21}^{+}$ by multiscale bootstrap (MB). The uniform density functions are given by the red lines.}
%\caption{Top row: Histograms of $p$-values of $H_{21}^{+}$ by ordinary bootstrap (BP) in the two variable cases. Bottom row:  Histograms of $p$-values of $H_{21}^{+}$ by multiscale bootstrap (MB) in the two variable cases. The uniform density functions are given by the red lines.}
\label{fig:twovar}
\end{figure}

\begin{figure}[!tb]
\begin{center}
\includegraphics[width=3in]{only6ss2.eps}
\caption{Scatterplots of empirical rejection probabilities of $H_{32}^{+}$ by ordinary bootstrap and those by multiscale bootstrap versus significance levels in the six variable cases. 
}
\end{center}
\label{fig:plots}
\end{figure}

The histograms of $p$-values of $H_{21}^{+}$ computed by ordinary bootstrap and those by multiscale bootstrap in the two variable cases are shown in Fig.~\ref{fig:twovar}. Similar histograms were obtained for the other conditions. 
Each of the histograms of $p$-values computed by multiscale bootstrap looked closer to the uniform distribution than by ordinary bootstrap. 
This implied that multiscale bootstrap provided better unbiased $p$-values. 

 In Fig.~3, we also show a scatterplot of empirical rejection probabilities by ordinary bootstrap  {\rm Prob}\{$p^{BP}(H_{32}$)$<$$\alpha$\} and those  by multiscale bootstrap {\rm Prob}\{$p^{MB}_3(H_{32}$)$<$$\alpha$\} versus significance levels $\alpha$ in the six variable cases. 
Plots for unbiased tests should be on the diagonal. 
That is, their rejection probabilities should be equal to the corresponding significance levels.
Most of the plots for ordinary bootstrap are far above the diagonal, indicating that ordinary bootstrap gave rather biased $p$-values and tended to reject reasonable hypotheses much more often than the nominal frequencies or significance levels.
In contrast, the plots for multiscale bootstrap are much closer to the diagonal. This showed that multiscale bootstrap provided much better unbiased $p$-values. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}\label{sec:conc}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We proposed a new procedure to evaluate statistical reliability of LiNGAM. 
Our procedure gives $p$-values of variable orderings estimated by LiNGAM and tells which orderings are more reliable. 
The utility of our procedure was demonstrated in the simulations. 
Future work would investigate how sensitive to non-smoothness of the boundaries of hypothesis regions our method is and how it is alleviated, although the simulations implied that it might be not very problematic. 
%The $R$-codes to implement our multiscale bootstrap procedure will be made available on the web. 

\subsubsection*{Acknowledgments}
This research was partially supported by JSPS Global COE program `Computationism as a Foundation for the Sciences'. 

% References
\bibliography{Shohei_Tex_Ref}
\bibliographystyle{splncs}

%\appendix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Proof of lemma \ref{lemma1}}\label{sec:proof1}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\end{document}













MBlingam_arXiv.bbl
\begin{thebibliography}{10}

\bibitem{Bollen89book}
Bollen, K.A.:
\newblock Structural Equations with Latent Variables.
\newblock John Wiley \& Sons (1989)

\bibitem{Pearl00book}
Pearl, J.:
\newblock Causality: Models, Reasoning, and Inference.
\newblock Cambridge University Press (2000)

\bibitem{Spirtes93book}
Spirtes, P., Glymour, C., Scheines, R.:
\newblock Causation, Prediction, and Search.
\newblock Springer Verlag (1993) (2nd ed. MIT Press 2000).

\bibitem{Shimizu06JMLR}
Shimizu, S., Hoyer, P.O., Hyv{\"a}rinen, A., Kerminen, A.:
\newblock A linear non-gaussian acyclic model for causal discovery.
\newblock J. Machine Learning Research \textbf{7} (2006)  2003--2030

\bibitem{Lehmann08Test}
Lehmann, E., Romano, J.:
\newblock Testing Statistical Hypotheses (3rd edition).
\newblock Springer (2008)

\bibitem{Efron93book}
Efron, B., Tibshirani, R.:
\newblock An Introduction to the {Bootstrap}.
\newblock Chapman \& Hall, New York (1993)

\bibitem{Felsenstein85Evol}
Felsenstein, J.:
\newblock Confidence limits on phylogenies: an approach using the bootstrap.
\newblock Evolution \textbf{39} (1985)  783--791

\bibitem{Fried99UAI}
Friedman, N., Goldszmidt, M., Wyner, A.:
\newblock Data analysis with {Bayesian} networks: A bootstrap approach.
\newblock In: Proc. Conf. on Uncertainty in Artificial Intelligence (UAI1999).
  (1999)  196--205

\bibitem{Efron96boot}
Efron, B., Halloran, E., Holmes, S.:
\newblock Bootstrap confidence levels for phylogenetic trees.
\newblock In: Proc. Natl. Acad. Sci. USA. (1996)  13429--13434

\bibitem{Hillis93SB}
Hillis, D., J.Bull:
\newblock An empirical test of bootstrapping as a method for assessing
  confidence in phylogenetic analysis.
\newblock Syst. Biol. \textbf{42} (1993)  182--192

\bibitem{Hall92book}
Hall, P.:
\newblock The bootstrap and Edgeworth expansion.
\newblock Springer-Verlag, New York (1992)

\bibitem{Shimo02SB}
Shimodaira, H.:
\newblock An approximately unbiased test of phylogenetic tree selection.
\newblock Systematic Biology \textbf{51} (2002)  492--508

\bibitem{Shimo08JSPI}
Shimodaira, H.:
\newblock Testing regions with nonsmooth boundaries via multiscale bootstrap.
\newblock J. Statistical Planning and Inference \textbf{138} (2008)  1227--1241

\bibitem{Hoyer07IJAR}
Hoyer, P.O., Shimizu, S., Kerminen, A., Palviainen, M.:
\newblock Estimation of causal effects using linear non-gaussian causal models
  with hidden variables.
\newblock Int. J. Approximate Reasoning \textbf{49} (2008)  362--378

\bibitem{Hyva01book}
Hyv{\"a}rinen, A., Karhunen, J., Oja, E.:
\newblock Independent component analysis.
\newblock Wiley, New York (2001)

\bibitem{Akaike74AIC}
Akaike, H.:
\newblock A new look at the statistical model identification.
\newblock IEEE Trans. Automat. Control \textbf{19} (1974)  716--723

\bibitem{Silva06JMLR}
Silva, R., Scheines, R., Glymour, C., Spirtes, P.:
\newblock Learning the structure of linear latent variable models.
\newblock J. Machine Learning Research \textbf{7} (2006)  191--246

\end{thebibliography}
