\documentclass[11pt]{article}
\usepackage{a4wide}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{natbib}
\usepackage{graphicx}
\usepackage{subfigure}
\usepackage{epsfig}\pssilent



\bibpunct{(}{)}{;}{a}{}{,}

\renewcommand\baselinestretch{1.4}



\begin{document}

\pagestyle{plain}
\setcounter{page}{1}
\setcounter{tocdepth}{2}


\begin{center}\Large\bf Supplementary Material to:

\begin{flushleft}
{\normalsize Exploration of Empirical-Bayes Hierarchical Modelling for 
the
Analysis of Genome-Wide Association Study Data}
\end{flushleft}

\normalsize \mdseries 
\begin{flushleft} $\mbox{Elizabeth A. Heron}^{\dagger  } 
\footnote{
Corresponding author. E-mail:
eaheron@tcd.ie\\
 \hspace*{4cm}Tel: +353(1)8962424}
\quad \mbox{Colm O'Dushlaine}^{\dagger}
\quad \mbox{Ricardo Segurado}^{\dagger}
\quad \mbox{Louise Gallagher}^{\dagger}
\quad \mbox{Michael Gill}^{\ddagger}$
\end{flushleft}

$\begin{array}{rl} \dagger  &\mbox{Neuropsychiatric Genetics
Research Group,}\\
&\mbox{Dept. of Psychiatry, Trinity College Dublin,}\\
&\mbox{Trinity Centre for Health Sciences, James's St.,}\\
&\mbox{Dublin 8, Ireland.}\\
\end{array}$


\end{center}



\vspace{0.5cm}


\section{Mixture Model and Bayesian Model Interpretations}

\label{bayes-freq}
Combining Equations 2.1 and 2.2 (Main Text), we
can consider the empirical-Bayes hierarchical model as a mixed model 
with both fixed and random
effects in the frequentist setting as given by:
\begin{eqnarray*}
\mbox{logit}(p_{mi}) = \alpha_{m} + Z_{m}\gamma X_{mi} +
X_{mi}\tau^{2}t_{mm},
\end{eqnarray*}
where \(Z_{m}\gamma X_{mi}\) would be considered a fixed effect and
\(X_{mi}\tau^{2}t_{mm}\) a random effect.
Alternatively, this model may be viewed in the Bayesian framework, where
the posterior distribution for the unknown parameters given the data is
given by:
\begin{eqnarray*}
P(\{\alpha_{m}\}, \{\beta_{m}\}, \{\gamma_{k}\}, \tau| \{y_{i}\})
&\propto &
\prod_{m=1}^{M}\prod_{i = 1}^{N}[p_{mi}^{y_{i}}(1 - p_{mi})^{1 -
y_{i}}]\\
&& \times \prod_{m=1}^{M}\prod_{k=1}^{K}
\left [\frac{1}{\sqrt{2\pi \tau^{2}t_{mm}}}
\exp{\left(
\frac{-(\beta_{m} - Z_{mk}\gamma_{k})^{2}}{2\tau^{2}t_{mm}}
\right)}\right]\\
&& \times
\pi(\tau^{2})
\prod_{m=1}^{M}\pi(\alpha_{m})
\prod_{k=1}^{K}\pi(\gamma_{k}).
\end{eqnarray*}
The term
\(\prod_{m=1}^{M}\prod_{k=1}^{K}
\left [\sqrt{2\pi \tau t_{mm}}^{1/2}
\exp{
\left( -(\beta_{m} - Z_{mk}\gamma_{k})^{2}(2\tau^{2}t_{mm})^{-1}\right)
}\right]\), which is the second-level regression (see Equation 2.2 (Main 
Text))
plays the role
of an informative prior distribution on the unknown effect size
parameters \(\beta\).
The term \(\prod_{m=1}^{M}\prod_{i = 1}^{N}[p_{mi}^{y_{i}}(1 -
p_{mi})^{1 -
y_{i}}]\) is the likelihood with the remaining terms, \(\pi(.)\)'s,
being prior distributions on the other unknown parameters. In this
context \(\tau^{2}\) and \(\{\gamma_{k}\}\) also have prior
distributions as they are also unknown parameters. A fully Bayesian
approach would allow for posterior distributional estimates of all the
unknown parameters of the model. For large numbers of markers,
estimating the unknown parameters of the model would not be feasible due
to the high dimensionality and complexity of the resulting model and
hence empirical-Bayes approaches have been developed. Presenting the
full Bayesian model here helps to clarify the role played by the various
components of the model.


\section{Empirical-Bayes Shrinkage Factors}

Starting with \(\gamma\), a weighted least squares estimator is given by:
\begin{eqnarray}
\label{gammas}
\hat{\gamma} = (Z^{T}DZ)^{-1}Z^{T}D\hat{\beta},
\end{eqnarray}
where \(D\) is the diagonal matrix of weights \(D = \mbox{Diag}(w_{1}, 
\ldots, w_{M})\), \(w_{m} = (V_{m} + \hat{\tau}^{2})^{-1}\), \(V_{m} = 
\mbox{SE}(\hat{\beta}_{m})^{2}\).
An approximately unbiased estimator of \(\tau^{2}\), the second-level 
variance, is given by:
\begin{eqnarray}
\label{A}
\hat{\tau}^{2} = 
\frac
{\sum_{m=1}^{M}w_{m}
\{
(M/(M-K))(\hat{\beta}_{m} - Z_{m}^{T}\hat{\gamma})^{2} - V_{m}
\}}
{\sum_{m=1}^{M}w_{m}}.
\end{eqnarray}
The process continues by iterating over Equations \ref{gammas} and 
\ref{A} until convergence is reached. An arbitrary initial guess is 
chosen for \(\hat{\tau}^{2}\). If \(\hat{\tau}^{2}\) becomes negative 
over the course of the iterations it is forced to be non-negative 
\((\hat{\tau}^{2} \ge 0)\). The shrinkage factor for each marker is then 
given by:
\begin{eqnarray*}
\hat{B}_{m} = ((M-K-2)/(M-K))V_{m}w_{m}.
\end{eqnarray*}

For the approximate empirical-Bayes confidence intervals the 
square of the standard 
error is given by (\citet{Morris1983}):
\begin{eqnarray}
s_{m}^{2} = V_{m}
[1-((M - \hat{k}_{m})/M)\hat{B}_{m}]
+ v_{m}(\hat{\beta}_{m} - Z_{m}^{T}\hat{\gamma})^{2},
\end{eqnarray}
where
\(\hat{k}_{m} = Mw_{m}[Z(Z^{T}DZ)^{-1}Z^{T}]_{mm}\) and \(v_{m} = 
(2/(M-K-2))\hat{B}_{m}^{2}(\bar{V} + \hat{\tau}^{2})/(V_{m} + 
\hat{\tau}^{2})\), which is an approximation for the variance of 
\(\hat{B}_{m}\) and the average variance is given by: \(\bar{V} = 
(\sum w_{m}V_{m})/ \sum w_{m}\). Thus a \(95\%\) approximate confidence 
interval for the empirical-Bayes estimate of the logarithm of the odds 
ratio is given by: 
\begin{eqnarray}
\hat{\beta}_{m; EB} \pm 1.96s_{m}.
\end{eqnarray} 

For further details see \citet{Morris1983}. 

\section{Simulation Strategy}

To simulate case-control data for the empirical-Bayes hierarchical 
model we begin by fixing the number of cases and controls to be 
\(N_{1}\) and \(N_{0}\),
respectively. Also the number of markers (SNPs), \(M\), is set. We then
generate data by starting at the top level of the hierarchical model
specified in Equation 2.2 (Main Text). The additional genotypic
covariate information summarised in the matrix \(Z\) is simulated as
either independent binary or multinomial variates. We note that in real
experimental covariate data, many of the entries in \(Z\) would be zero
and we try
to mimic this with the choice of parameters for either the Binomial or
Multinomial distributions from which we simulate.
We choose one of
the columns of \(Z\) to contain the most informative additional
covariate information and hence the most influential as regards a SNP
being causative or associated. This is a reasonable assumption as it
will often be the case that one or more of the additional types of
genotypic data will be more informative than the others (a rationale
used by \citet{ChenWitte2007} in the development of the
weighting function in their empirical estimators).

At this stage a number of SNPs (at least 1) are chosen to be associated,
by this we mean that in our simulation studies the odds ratios for these
particular SNPs will be large relative to all other SNPs simulated. The
coefficients \(\gamma\) of the second-level regression, given in
Equation 2.2 (Main Text), are fixed as constants, with the
coefficient of the covariate that was chosen to be most informative set 
at a
higher value with respect to all others. The mean of
the second-level regression is then given by \(Z \gamma \).
A variance-covariance matrix is constructed: \(\tau^{2}T\). This matrix
represents the residual variation that remains in the first-level
coefficients (the \(\beta_{m}\)'s) after the second-level covariate
information has been incorporated. For computational ease we do not
simulate any correlations between the SNPs and so \(T\) can be set
as the identity matrix. \(\tau^{2}\) is not chosen to vary across SNPs
but rather it is constant and thus we are, in this case, modelling equal
residual variability across the SNPs. Given the mean and variance we can
now simulate \(\beta_{m}\)'s from the appropriate Normal distributions.
This completes the first part (second-level) of the simulation.

The second part of the simulation study involves simulating realistic
genotypes for the \(N\) individuals for each of the \(M\) SNPs. We
start with the MAFs, which can be either
fixed at constant values, or they can be randomly drawn from a Uniform
distribution on \((0, 0.5)\). In order to simulate the genotypes it is
also necessary to consider the prevalence of the disease, which is
defined to be the probability of disease, and which we fix at a constant
value.
Simulation of the genotypes for each marker is considered for cases and
controls separately. Specifically, for each marker, for each individual,
we simulate the genotypes using a multinomial distribution with 
parameters
\(P(aa_{m}|\mbox{case})\), \(P(Aa_{m}|\mbox{case})\) and 
\(P(AA_{m}|\mbox{case})\),
if the
individual is
a case and \(P(aa_{m}|\mbox{control})\), \(P(Aa_{m}|\mbox{control})\) 
and
\(P(AA_{m}|\mbox{control})\),
if the individual is a control, for the three genotypes \(aa_{m},
Aa_{m}\) and \(AA_{m}\). In order to calculate these
six probabilities for each marker it is necessary to specify the genetic
model, e.g, dominant, additive, etc., and to calculate the baseline risk
of disease for each marker, \(\alpha_{m}\), as all other parameters have
been set.
Further details of how these are calculated are given below. Thus we 
now have a strategy for simulating a set of cases and controls with each 
of their genotypes for the \(M\) markers. All data simulation and 
analysis is carried out in the statistical computing language: 
\citet{R}.

\subsection{Probability of an Individual Being a Case Given Their Genotype}

The probability that an individual has the genotype \(aa_{m}\) at marker 
\(m\) given that they are a case is given by the following:
\begin{eqnarray*}
P(aa_{m}|\mbox{case}) 
& = & \frac{P(\mbox{case}|aa_{m})P(aa_{m})}{P(\mbox{case})}\\
& = & \frac{P(\mbox{case}|aa_{m})(1-f_{m})^{2}}{p},
\end{eqnarray*}
where \(f_{m}\) is the MAF (allele \(A\) is the minor allele) at the 
marker 
\(m\) and \(p\) 
is the prevalence of the disease. Both of these are fixed in the 
simulation study. \(P(\mbox{case}|aa_{m}) = \exp{(\alpha_{m})}/(1 + 
\exp{(\alpha_{m})})\), and it is necessary to solve for \(\alpha_{m}\), 
the baseline risk of disease for each marker \(m\). Similarly, 
expressions for the probabilities of the other genotypes given 
case-control status can also be derived. Solving for \(\alpha_{m}\) 
allows for all of these to be calculated. 

\label{appendix probabilities}

\subsection{Baseline Risk for the Dominant Model}

In order to obtain an expression for the baseline risk of disease for 
each marker, \(\alpha_{m}\), the genetic model must be taken into 
account. For a dominant model, using the logistic regression given in 
Equation 2.1 (Main Text) and assuming Hardy-Weinberg 
Equilibrium, the 
prevalence of disease is given by:
\begin{eqnarray}
p = 
\frac{e^{\alpha_{m}}}{1 + e^{\alpha_{m}}}
(1 - f_{m})^{2}
+
\frac{e^{\alpha_{m} + \beta_{m}}}{1 + e^{\alpha_{m} + 
\beta_{m}}}(2f_{m} - f_{m}^{2}),
\label{prob dominant}
\end{eqnarray}
To obtain an expression for \(\alpha_{m}\) we re-write Equation 
\ref{prob dominant} as a quadratic equation in \(e^{-\alpha_m}\). To 
simplify notation we let \(x = e^{-\alpha_m}\), \(y = e^{-\beta_{m}}\), 
and \(f = f_{m}\). The quadratic equation is then given by:
\begin{eqnarray}
x^{2}py + xB + p-1 = 0,
\label{quadratic dominant}
\end{eqnarray}
where \(B = (2f - f^{2} - 1)y + (f^{2} - 2f) + (p + yp)\). Solving 
Equation \ref{quadratic dominant} and considering only the positive 
root, results in 
\begin{eqnarray}
\alpha_{m} = 
-\log{\left(
\frac
{(-B + \sqrt{(B^{2} + 4(1-p)yp)})}
{2yp}
\right)} 
\label{alpha dominant}
\end{eqnarray}
as an expression for the baseline risk of disease for the dominant 
model.

\subsection{Baseline Risk for the Additive Model}
For the additive model, it is necessary to solve a cubic equation given 
by:
\begin{eqnarray*}
&& x^{3}y^{3}p - x^{2}[y^{3}(1-f^{2}-p) + y^{2}(2f-2f^{2}-p) + 
y(f^{2}-p)] 
\\
&& - x[(y^{2} + y)(1 - f^{2}) + 2(y^{2} + 1)(f - f^{2}) + (y+1)f^{2} - 
p(y^{2} + y + 1)] 
\\
&& - (1 + 2(f - f^{2}) -p) -p = 0,
\end{eqnarray*}
where \(x, y, f, p\) are as in Equation \ref{alpha dominant}.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thebibliography}{99}

\bibitem[Chen and Witte(2007)]{ChenWitte2007}
         Chen, G.K., and Witte, J.S. (2007),
         ``Enriching the Anlaysis of Genomewide Association Studies with
         Hierarchical Modeling," {\it The American Journal of Human
         Genetics}, 81, 397-404.

\bibitem[Morris(1983)]{Morris1983}
         Morris C.N. (1983),
         ``Parametric Empirical-Bayes Inference: Theory and 
           Applications," {\it
           Journal of the American Statistical Association}, 78(381), 
	 47-55.

\bibitem[R(2009)]{R}
         R Development Core Team (2009),
         ``R: A Language and Environment for Statistical Computing,"
          {\it R Foundation for Statistical Computing}, Vienna,
          Austria, ISBN: 3-900051-07-0,
          http://www.R-project.org.


\end{thebibliography}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage



\newpage

\begin{figure}[!ht]
\begin{center}
        \subfigure[]{
          \includegraphics[height = 6cm, width = 7.2cm, angle =
           0]{Heron_SM_Fig1a}
        }
        \subfigure[]{
        \includegraphics[height = 6cm, width = 7.2cm, angle =
         0]{Heron_SM_Fig1b}
        }
\end{center}
\caption[MSE plots]{
(a) MSE for both the MLEs and the empirical-Bayes estimators (EB-HM)
across a range of 100 random simulations, when the same \(Z\) matrix is
used
in the simulation and analysis phases.
Also plotted is the proportion of MAFs \(< 5\%\)  in each simulation.
(b) MSE for both the MLEs and the empirical-Bayes estimators
across a range of 100 random simulations, when the alternative reduced,
noisy \(Z\) matrix is used in the analysis phase.
Also plotted is the proportion of MAFs \(< 5\%\)  in each simulation.
}
\label{MSE Plots}
\end{figure}





\newpage

\begin{figure}[!ht]
\begin{center}
        \subfigure[]{
          \includegraphics[height = 6cm, width = 7.2cm, angle =
           0]{Heron_SM_Fig2a}
        }
        \subfigure[]{
        \includegraphics[height = 6cm, width = 7.2cm, angle =
         0]{Heron_SM_Fig2b}
        }
          \ \\
        \subfigure[]{
          \includegraphics[height = 6cm, width = 7.4cm, angle =
           0]{Heron_SM_Fig2c}
        }
        \subfigure[]{
        \includegraphics[height = 6cm, width = 7.4cm, angle =
         0]{Heron_SM_Fig2d}
        }
\end{center}
\caption[Top 200 plots]{
(b) The top 200 markers with the largest EB-HM odds
ratios together with their approximate 95\% confidence
intervals for (a) CAD, (b) CD, (c) RA, and (d) T2D.

}
\label{Top 200 Plots}
\end{figure}

\newpage








\end{document}

