

\documentclass[12pt]{article}
\title{Notes for Different Graphical Gaussian Models}
\author{Guiyuan Lei\footnote{Guiyuan.Lei@ncl.ac.uk}}
\date{September, 2006}
\begin{document}
\maketitle
\section{HdBCS: Bayesian Covariance Selection in High-Dimensions}
The model of Bayesian Covariance Selection in High-Dimensions 
\cite{Dobra_2004} is designed to perform covariance selection for datasets with tens or possibly hundreds of thousands of variables. 
\subsection{The algorithm}
The main idea of HdBCS is to identify good regression models for each variable and to combine these models in a joint multivariate normal distribution using the chain rule (which is essentially equivalent to constructiong a DAG or Bayesian network). The algorithm starts with dependency network which is built from sets of regressions seperately, then generates an appropriate ordering of the variables to underlie a compositional network. With this ordering, variable selection for each linear model for the ordered variables produces a DAG. Finally, by removing arrows and marrying all parents of any node in this DAG, an undirected graph is delivered.\\

For analysis and selection of regressions, the model utilizes the inverse Wishart prior, applies for forward/backward varialbe selection with respect to the defined set of candidate predictors, and score each regression by the posterior model probabilities. See section 3 in the paper \cite{Dobra_2004} for detail. \\
%equation (4) is marginal model likelihood

Once the predictors of each variable are obtained, the order of all variables should be generated to construct compositional network. The strategy is to: the ordering aims to maximize the resulting posterior probability on the resulting model in an iterative process that ,at each step, decreases the overall score (equation 5 in the paper) the least. See section 4.1 for fully described construction.

\subsection{The software}
There is one software implemented for this model. Codes can be downloaded from 
http://www.stat.duke.edu/$\sim$adobra/hdbcs.html\\
%http://www.stat.duke.edu/~adobra/hdbcs.html
%HdBCS: Bayesian Covariance Selection in High-Dimensions
%Adrian Dobra (adobra@stat.duke.edu)

The search procedure has three separate steps: filtering, generating good starting models and improving the starting models until convergence. \\ 

The fisrt step is to select a relatively rich set of possible predictors for each variable. The search strategy is forward/backward procedure. This procedure is for each variable seperately. Starting with random selected predictors for each variable $x_i$, consider the edges between all other variables $x_j, j=[-i]$ with this variable, if $x_j$ is already predictor of variable $x_i$, propose to delete $x_j$ as predictor of $x_i$, otherwise, propose to add $x_j$ as its predictor. Calculate the ratio of posterior probability between proposed model and current model (one predictor difference between these two models), and accept/reject the proposed model according to the accept propobility. Repeat the above process for many iterations to get set of possible predictors for each variable. The pseudo code for the search procedure can be describled as following:\\

%see gibbsreg.cpp
For each variable $x_i$, 
\begin{enumerate}
\item Initial model, randomly select predictors of $x_i$ from all other variables $x_j, j=$\{$1,\cdots,p$\}$\setminus$\{$i$\} 
\item Iteration step: repeat the following process for a fixed times:
 \begin{enumerate}
 \item[] For all other variables $x_j, j=$\{$1,\cdots,p$\}$\setminus$\{$i$\}, randomly select $x_j$ (each variable can be selected once)
       \begin{enumerate}
        \item Propose new model from current model    
          \begin{enumerate}
             \item[] If ($x_j$ is not currently predictor of $x_i$), calculate the log posterior probability with $x_j$ is   predictor of $x_i$
             \item[] else, calculate the log posterior probability calculate the log posterior probability with $x_j$ is not predictor of $x_i$
           \end{enumerate}
         \item Calculate accept propobility and accept/reject proposed model, if $x_j$ is a predictor in new model, update the frequency of $x_j$ as predictor of $x_i$ % varFreq[Predictors[i]]+=1
        \end{enumerate}
 \end{enumerate} %end of "repeat the following process"
\item For all other variable $x_j, j=$\{$1,\cdots,p$\}$\setminus$\{$i$\}, if the frequency $x_j$ as predictor of $x_i$ is greater than threshold, set $x_j$ as predictor of $x_i$ in filtered model (one dependency network) %varFreq[]
\end{enumerate}
%=======
%in gibbsreg.cpp
%void gibbsSamplingReg(int targetGene,int* varFreq)
%if(1==VarInModel[Predictors[i]])
%	    {
%               varFreq[Predictors[i]] += 1;
%	    }
%========

In second step, first generate dependency network which is based the filtered model in step1, the possible predictors for initial model are selected from the final model in step1 (in step1, all other variables [$-i$]=\{$1,\cdots,p$\}$\setminus$\{$i$ \} are possible predictors of $x_i$), in addition, the best model is recorded (in step1, final model is average model); then DAG is generated from the best model. In this step, generate several DAGs by running several times of above procedure. For generating compositional network (DAG) from dependency network, see section 4.1.2 in the paper \cite{Dobra_2004}.\\
%==========
%in startpoint.cpp
%//this graph is for keeping the set of possible predictors 
%LPList predictorGraph = InitGraph(NumberOfGenes);
%ReadGraph(predictorGraph,infile);

%in variousfunctions.cpp
%int gibbsSamplingReg(int targetGene,LPList predictorGraph,int* bestModel)
%bestModel[Predictors[i]] = VarInModel[Predictors[i]];
%==========

In third step, the starting points/models (which are generated in step2) are improved in a procedure (simulated annealing) that converges to local models of the posterior distribution over the space of DAGs. \\

Given possible set of predictors $S_i$ for each variable $x_i$ and DAG generated in step2. With order of variables for given DAG, switch variable $x_i$ in position $equation$ with variable $x_j$ in position $equation+1$, reconsider the regression equations of these two varialbes. The search for regression equation of variable $x_i$ is similar as that of in step1, but the predictors are restricted to the intersection of set $S_i$ and successor of $x_i$. 
%$S_i$ is recorded in list predictorGraph which is read from file *.predictors, in slave() sub-program. predictorGraph is used in sub-program gibbsSamplingReg() to select predictor once and only once (see ChooseNextGene()), it is important that the vector mymodel is reset to zero after each iteration in gibbsSamplingReg().
%
Calculate the posterior probabilies for variable $x_i$ and $x_j$, also calculate the difference $\delta$ between these posterior probabilies and that of former regression equations for variable $x_i$ and $x_j$. Calculate $\delta$ for each varialbe $x_i, i=0,1,\cdots,NumberOfGenes-1$ anb sort these $\delta$ in a list $eq$. \\
%The above process is ``generate candidates for ALL regression equations`` in masterGraph(), where the search process for regression equation is the case of "`"FIND_CANDIDATE' in slave()

Next, select which equation (and its adjacent equation+1) to be switched order. The weights $w[i]$ of each pair of order switching are normalized as following: 
\begin{enumerate}
             \item $w[i]=w[i]-w[nmax-1]$, where $nmax$ is number of pair switching, the weights are in increasing order
             \item $w[i]=pow(exp(w[i]),anneal)$, where $anneal$ is annealing parameter, the minimal value is $1$, maximum valuse is set to be $3$
             \item $cumw[i]=\sum_{j=0}^{i-1}x[j], i=1,2,\cdots,nmax$
             \item $cumw[i]=cumw[i]/cumw[nmax]$
\end{enumerate}
The pseudo code for selecting equation by weights $cumw[i]$ is as following:
\begin{enumerate}
             \item Initialize: $q1=0;q2=nmax;$ and generate a random number $s$
             \item While($q1+1!=q2$) 
                     \begin{enumerate}
                         \item[]$ q=(q1+q2)/2$;
                         \item[] If($cumw[q]<s$), $q1=q$;
                         \item[] else, $q2=q$      
                     \end{enumerate}
\end{enumerate}
The returned $q1$ is the selected equation, denote it as $chosenEquation$ (this is the position of selected equation/variable in the order)
%SelectEquation()
%NormalizeWeights() and WeightedSampling()
and update order and regression equations for $x_i$ and $x_j$ in position $chosenEquation$ and $chosenEquation+1$ in the order.\\
%case ``UPDATE_ORDERING`` and ReplaceRegressions()

Then, for generate next DAG, only need to calculate $\delta$ for three pairs of ordering switch, that is $chosenEquation-1$ (and its adjacent position $chosenEquation$), $chosenEquation$ (and its adjacent position $chosenEquation+1$) and $chosenEquation+1$ (and its adjacent position $chosenEquation+2$).
%For graphnumber from $1$ to $nGraphsToGenerate$, the search for regression equation for two switched variable is done seperately (for graphnumber 0, the case ``FIND_CANDIDATE``, 	the search for regression equations for two switched variable is done together), that is to execute ``FIND_CANDIDATE_UPDATE`` twice.
After updating the $\delta$ for new order, select equation to switched with variable in its adjacent position, update order and regression equation for next new DAG.\\
%UpdateEquationsScores(), SelectEquation(), ``UPDATE_ORDERING``, ReplaceRegressions()

Repeat the above procedure for $nGraphsToGenerate$ times and record these DAGs/models every $whenToSave$ iterations.\\

After all these three steps, an average network is generated from all statistically significant models during the simulated annealing procedure. Find the model with best posterior probability $w$, all other recorded (for example, recorde every 1000 iterations) models with porstiror probability $w_i$
\begin{equation}
\exp{(w_i-w)}>\epsilon
\end{equation}
will be considered in averaging network. When add one graph to averaged graph, the weight of each edge $w_i$(if this edge exists in current graph) is added as $b+\log{(1+\exp{(a-b)})}$ (assuming $a<b$, otherwise just exchange the value of $a$ and $b$) if the weight of average edge is $a$ (initialized as $0$) and the weight of added graph is $b$. Note that $a$ and $b$ are log posterior probability and $b+\log{(1+\exp{(a-b)})}=\log{(\exp{(a)}+\exp{(b)})}$. Meanwhile, the weight of average graph $GrandTotal$ is added as same formula $b+\log(1+\exp(a-b))$ but the value of $a$ is set to be very small number, like $-99999999999999999.0$. At last, when all statistically significant graphs are considered, the weight of each edge is $\exp{(w_i-GrandTotal)}$.

%MatrixAddNumbers() in sparsematrix.cpp, If a<b, ruturn b+log(1+exp(a-b));

%averagedags.cpp, in main()
%double GrandTotal = -99999999999999999.0;
%GrandTotal = MatrixAddNumbers(GrandTotal,DagWeight[i]);

%in SaveMatrix(LPMatrix matrix,LPMatrix isBest,double GrandTotal)
%exp(p->value-GrandTotal)) for the average weight of each edge


%\addcontentsline{toc}{chapter}{\numberline{}Bibliography}
%\bibliography{GGM}
\begin{thebibliography}{}

  \bibitem{Dobra_2004}
     Adrian Dobran, Chris Hans, Beatrix Jones, Joseph R. Nevins, Guang Yao \& Mike Wes. 
     {Sparse graphical models for exploring gene expression data}
     \textit{Journal of Multivariate Analysis} \textbf{90} (2004), 196--212.
\end{thebibliography}
\end{document}


