\documentclass[11pt]{article}
\usepackage{amsmath,amsfonts,Sweave,natbib}
\usepackage{hyperref}
 \usepackage[utf8]{inputenc}
\usepackage[margin=1in]{geometry}



\newcommand{\proglang}{}
\newcommand{\pkg}[1]{ {\bf #1}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% declarations for jss.cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% almost as usual
%\author{Arend Voorman\\University of Washington 
 %     others}
\title{skatMeta: an R Package for meta analyzing region-based tests of
rare DNA variants}

%% for pretty printing and a nice hypersummary also set:
%\Plainauthor{Arend Voorman, more people} %% comma-separated
%\Plaintitle{skatMeta: an R Package for meta analyzing region based tests} %% without formatting
%\Shorttitle{skatMeta} %% a short title (if necessary)

%% an abstract and keywords

%\Keywords{meta-analysis, genetics, \proglang{R}}
%\Plainkeywords{meta-analysis, genetics, R} %% without formatting
%% at least one keyword must be supplied
%% publication information
%% NOTE: Typically, this can be left commented and will be filled out by the technical editor
%% \Volume{50}
%% \Issue{9}
%% \Month{June}
%% \Year{2012}
%% \Submitdate{2012-06-04}
%% \Acceptdate{2012-06-04}

%% The address of (at least) one author should be given
%% in the following format:
%\Address{
%  Arend Voorman\\
%%  Department of Biostatistics\\
% University of Washington\\
%  F-649, Health Sciences Building, Box 357232 \\
%  Seattle, WA 98195-7232
%  E-mail: \email{voorma@uw.edu}\\
%}
%% It is also possible to add a telephone and fax number
%% before the e-mail in the following format:
%% Telephone: +43/512/507-7103
%% Fax: +43/512/507-2851

%% for those who use Sweave please include the following line (with % symbols):
%% need no \usepackage{Sweave.sty}

%% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}
\maketitle
 %\VignetteIndexEntry{skatMeta}
 %\VignetteKeywords{skat,skatMeta}
\begin{abstract}
Region-based tests are becoming a popular tool for analyzing rare
genetic variants. In order for these tests to have adequate power, it is often necessary
to meta-analyze information from multiple contributing cohorts, where consent
restrictions make it difficult or impossible to share individual level data. We
present the \proglang{R} package \pkg{skatMeta} for meta-analyzing region based tests,
such as SKAT, SKAT-O, and burden tests, as well as single variant
tests. The package can accommodate continuous, binary, and survival
outcomes for unrelated individuals, and continuous outcomes for
related individuals. We also provide convenient functions for
conditional analyses and parallelization.
\end{abstract}


%\setkeys{Gin}{width=0.6\textwidth}
\section{Introduction}
In association studies of rare genetic variation it is common to pool
variants across a gene, region, or pathway to improve statistical
power and performance of tests. In order for these test to have
adequate power to detect rare variants with modest effects, it is often
necessary to use data from multiple contributing cohorts using a meta-analysis. While meta-analyzing single-variant
effect estimates and their standard errors is
relatively straightforward, implemented in e.g. METAL, meta-analyzing
region-based tests is somewhat more complicated.  One popular region
based test is SKAT, which is based on a weighted sum of single-variant score
tests, using the covariance matrix of the scores to
calculate $p$-values \citep{wu2011rare} . Meta-analyzing SKAT
requires meta-analysis of score vectors and their covariance
matrices. However, meta-analyzing covariance matrices is cumbersome,
as they are 2-dimensional structures of various sizes and cannot be
easily stored for a whole study in an
array. The purpose of the \pkg{skatMeta} package is to facilitate meta analysis of region-based score tests,
such as SKAT, through convenient data structures for handling the
components of the tests. There is some software available for similar
analyses, for instance the MASS program
(\url{http://dlin.web.unc.edu/software/mass/}). Our implementation is
written in \proglang{R} \citep{R} and supports continuous, binary and survival
outcomes for unrelated individuals, continuous outcomes for related
individuals, and performs individual variant tests, as well as
burden tests, SKAT and SKAT-O.

Use of the package proceeds in
two steps. The first is running the cohort level analyses, with either
\verb=skatCohort= for unrelated individuals, or \verb=skatFamCohort=
for family data. These functions need only be called once
in each cohort for all genes under consideration. The output is a an
\proglang{R} `list' object
where each element corresponds to a gene, and contains the scores and
MAFs for that gene, as well as  a matrix of the covariance between the scores at all
pairs of SNPs within the gene. This list of output is labeled as a
\verb=skatCohort= object, is
saved as an \verb=.Rdata= file, and passed to a central location. The
second step combines the output from the cohort-level analyses,
i.e. doing meta-analysis. The meta-analysis functions can create
several different tests including single variant, T1 count
\citep{li2008methods}, Madsen-Browning \citep{madsen2009groupwise}, and
of course SKAT  \citep{wu2011rare}.  

Meta-analysis functions can be run on a single cohort alone, and this
will produce cohort level test results. These can be useful for
quality control. But the meta-analysis across cohorts does require the full output from \verb=skatCohort=/\verb=skatFamCohort=; it cannot be done with just cohort level test results. 

The rest of this document is organized as follows. Section
\ref{sec:basic} details the functions used at the cohort level for
continuous outcomes in unrelated individuals
(Section~\ref{sec:unrelated}), and related individuals
(Section \ref{sec:family}). Section~\ref{sec:meta} describes the
functions used at the meta analysis level for SKAT (Section~\ref{sec:skat}),
burden tests (Section~\ref{sec:burden}) and single SNP tests
(Section~\ref{sec:ss}). Section~\ref{sec:other} describes
modifications for binary and survival
outcomes. Section~\ref{sec:conditional} describes how to perform
conditional analyses. Section \ref{sec:bigdata} gives some guidance for parallelization on
large datasets, and finally Section~\ref{sec:method} describes
the mathematical details of the implementation. All sections contain
example code which should be directly executable in \proglang{R} using the
package \pkg{skatMeta} available at \url{http://cran.r-project.org/web/packages/skatMeta/}

\section{Cohort level analysis}
\label{sec:basic}
We are interested in testing the association
of groups of SNPs with a particular outcome. A fundamental ingredient for region-based tests is a common
definition across cohorts of which variants are in which regions. For convenience, we will
refer to these regions of SNPs as genes, though in principal they can be any unit of
aggregation. All functions in the
package require a SNP Information file (argument \verb=SNPInfo=), which
links SNPs to genes. By default, the package looks within SNP
Information file (a data frame) for the fields
\verb=Name= and \verb=gene= for SNP and gene identification respectively. If other
fields are desired to indicate SNP names or genes, these must be
specified explicitly via the arguments \verb=snpNames= and \verb=aggregateBy=.

Using \verb=skatCohort=, users are strongly recommended to include all
SNPs in the analysis, even if they are monomorphic in one study. This is for two
reasons; firstly, monomorphic SNPs provide information about MAF
across all studies; without providing the information we are unable to
tell if a missing SNP data was monomorphic in a cohort, or simply
failed to genotype adequately in that cohort. Second, even if some
SNPs will be filtered out of a particular analysis (e.g., because
they are intronic or common) constructing skatCohort objects
using all available SNPs allows greater flexibility in the meta-analysis.

\subsection{Unrelated subjects}\label{sec:unrelated}
To create a \verb=skatCohort= object containing one cohort's contributions, use e.g.
\begin{verbatim}
skatCohort(Z=Z, formula =nullModel, 
          SNPInfo=SNPInfo, data=pheno)
\end{verbatim}
where;
\begin{itemize}
\item \verb=Z= is a NSNP rows $\times$ NSAMPLE columns matrix
  containing the additively coded (typically 0-1-2 for number of minor
  alleles) genotypes, with column names which match the SNP names used
  in the SNPInfo file. Rows do not have to be named, but must be
  in the same order as in the phenotype file (the \verb=data=
  argument).  The \verb=Z= matrix
  can contain missing genotypes (coded as NA); missing genotypes are
  imputed to the mean value of the genotype in that cohort. Because
  there is no allele-checking functionality, all genotypes should be
  coded the same way for all cohorts before running. 

\item \verb=formula= is a formula object for the null model of the form
\verb=outcome ~ covariates=. In the example below, we use the formula \verb=y~sex+bmi=
in the example below to analyze an outcome \verb=y=, adjusted for covariates \verb=sex= and \verb=bmi=.

\item \verb=SNPInfo= is a matrix or data-frame with columns titled
  \verb=Name= and \verb=gene=. These columns specify the SNPs, and
  their associated genes, to be used in the analysis. One can use
  different columns to designate SNP names and genes with the arguments \verb=Name=
and \verb=aggregateBy= respectively.

\item \verb=data= is a data frame containing the phenotype information. Phenotypes and genotypes must have the same number of rows and be ordered in the same way, because \verb=skatCohort= matches them simply on row order and not by any form of ID variable. The phenos data frame must contain no missing data in the columns used in the analysis.  
\end{itemize}

For example, analysis at a single cohort, using the example data in
the package, might proceed as below.
<<>>=
rm(list=ls())
library(skatMeta)
######### load example data:
# contains SNPInfo, phenotyes (pheno1, pheno2)
# genotypes (Z1,Z2), and pedigree information (kins) for pheno2
data(skatExample)
ls()

#Perform cohort-level analysis:
c1 <- skatCohort(Z1, y~bmi+sex, SNPInfo = SNPInfo, data = pheno1)

###save the output, which can be passed to a central location.
cohort1.out.file <- tempfile()
save(c1, file = cohort1.out.file)
@ 

In this section we described analysis for continuous outcomes.  Binary
and survival outcomes for unrelated individuals require only slight
modification, and are described in Section~\ref{sec:other}.

\subsection{Related subjects}\label{sec:family}
This analysis proceeds in essentially the same was as for related
subjects, but uses \verb=skatFamCohort= instead of
\verb=skatCohort=. In addition to the arguments required by
\verb=skatCohort=, \verb=skatFamCohort= requires pair-wise kinship
information in matrix format.

As in \verb=skatCohort=, in \verb=skatFamCohort= phenotypes and
genotypes must have the same number of rows and these rows must be in the same order.  The function is matching them based on row order and not by any type of IDs.

<<>>=
c2 <- skatFamCohort(Z2, y~bmi+sex, SNPInfo = SNPInfo, fullkins = kins, 
                    data = pheno2)
cohort2.out.file <- tempfile()
save(c2, file = cohort2.out.file)
@ 

Family structure information is usually stored on disk in a `makeped'
linkage format. To convert this format into to the required kinship matrix, the
\verb=makekinship= function in the \pkg{kinship2} \proglang{R}!
package is recommended; see these packages documentation for
details \citep{kinship2}.  
 
\section{Meta analysis}\label{sec:meta}
The results of the one or more cohort level analyses can be
meta-analyzed using the functions \verb=skatMeta=, \verb=skatOMeta=, \verb=burdenMeta=, or
\verb=singlesnpMeta= for SKAT, SKAT-O, burden tests and
individual SNPs respectively. Regardless of the individual
analysis (related or unrelated individuals with continuous, binary or
survival outcomes) the meta analysis procedure is the same. Each of
these meta-analysis functions has a similar
syntax, and contains the arguments

\begin{itemize}
\item \verb=...=  One or more \verb=skatCohort= objects to be meta-analyzed.
  \item\verb=SNPInfo= The SNP Info file, listing genes and SNPs to be
    included in the meta analysis. This is in the same format as the
    SNPInfo file in \verb=skatCohort=, and should use the same naming
    convention for SNPs and genes.
   \item\verb=snpNames= The field of SNPInfo where the SNP identifiers are found. Default is \verb=Name=
  \item\verb=aggregateBy= The field of SNPInfo on which the skat results were aggregated. Default is \verb=gene=. For single snps which are intended only for single variant analyses, it is recommended that they have a unique identifier in this field.
 \item\verb=mafRange= A range of minor allele frequencies to be included
     in the analysis.  These are based on the pooled MAFs over all
     cohorts. By default this is c(0,0.5), which includes all
     SNPs.  
    \item\verb=verbose= logical, whether or not to print progress
     bars. Defaults to FALSE.
     
\end{itemize}

Each of the functions permits additional arguments, which we detail in
their corresponding sections. 

\subsection{Meta analyzing SKAT}\label{sec:skat}
The SKAT test is a weighted sum of individual score statistics
$$Q = \sum_j w_jU_j^2$$
where $w_j$ is a weight and $U_j$ is the score statistic for for the
association between phenotype and variant $j$.
Details are given in
Section~\ref{sec:method}. In the function \verb=skatMeta!= these weights are
specified via the \verb=wts= argument, which
gives either a function to be applied to the minor allele frequencies
or a column of the SNP Information file which can be coerced to a
numeric vector.

The additional argument \verb=method= allows the user to specify the
method of $p$-values calculation. The default is
\verb=method=`saddlepoint'=, which appears to work well in
practice. See the documentation of \verb=pchisqsum= and
Section~\ref{sec:method} for more details. 

For example, below is a meta analysis of the SKAT test.  Here we have
use the default  testing weights \verb=wts = dbeta(maf,c(1,25))=,
often referred to as the `Wu' weights,
which are also the default in the \pkg{SKAT} package \citep{skat}.
<<>>=
load(cohort1.out.file)
load(cohort2.out.file)

meta.results <- skatMeta(c1, c2, SNPInfo = SNPInfo)

head(meta.results)
@ 
Here, \verb=skatMeta= returns the following information
\begin{itemize}
  \item{\verb=gene=} The gene name, or other unit of aggregation.
  \item{\verb=p=} The $p$-value from the SKAT test
  \item{\verb=Qmeta=} The SKAT $Q$ statistic, defined as $\sum_j w_j
   U_j^2$, where $w_j$ is the weight given to SNP $j$, and $U_j^2$
    is associated score statistic.
   \item{\verb=cmaf=} The cumulative minor allele frequency. That is
     $\sum_j \mathrm{MAF}_j$, where $\mathrm{MAF}_j$ is the minor
     allele frequency of variant $j$, and the sum is over all variants
     in the gene.
   \item{\verb=nmiss=} The number of missing SNPs.  For a gene with a
     single SNP this is the number of individuals which do not
     contribute to the analysis, due to cohorts that did not report
     results for that SNP. For a gene with multiple SNPs, \verb=nmiss=
     is summed over the gene.
   \item{\verb=nsnps=} The number of SNPs in the gene.
   \end{itemize}
      
 It will often be useful to examine results within cohort before meta
analysis.  This can be done by `meta-analyzing' a single cohort:

<<>>=
cohort1.results <- skatMeta(c1, SNPInfo = SNPInfo)
@

\subsection{Meta analyzing burden tests}\label{sec:burden}
 Another commonly used family of tests are the `burden'
tests which regress the phenotype on a weighted sum of genotypes, within each
gene.
The score test for a weighted sum of genotypes has the form 
$$T = \sum_j w_j U_j,$$ where $w_j$ is a weight for SNP $j$ and $U_j$
is the score for SNP $j$.  For instance if $MAF_j$ is the minor allele
frequency for SNP $j$, the Madsen-Browning \citep{madsen2009groupwise}
test uses $w_j = (MAF_j(1-MAF_j))^{-1}$, while the
and the T1 count test \citep{li2008methods} uses $w_j = \mathbf{1}(MAF_j < 0.01)$. 

These types of tests can be computed with the
\verb=burdenMeta= function. This has the argument \verb=wts=, which
like \verb=skatMeta=, gives either a function to be applied to the minor allele frequencies
or a column of the SNP Information file which can be coerced to a
numeric vector.

To perform the T1 test one could use the command
<<>>=
meta.t1.results <- burdenMeta(c1, c2, wts = function(maf){maf < 0.01}, 
                              SNPInfo = SNPInfo)
@ 
Equivalently, we could use constant weights and limit the analysis to those SNPs with $MAF < 0.01$:
<<>>=
meta.t1.results <- burdenMeta(c1, c2, wts =1, 
                              mafRange = c(0,0.01), SNPInfo = SNPInfo)
@ 
Likewise, the Madsen-Browning test could be performed with the command
<<>>=
meta.mb.results <- burdenMeta(c1, c2, wts = function(maf){1/(maf*(1-maf))}, 
                              SNPInfo = SNPInfo)
@ 
Regardless of the weights used, these tests give output of the form
<<>>=
format(head(meta.t1.results),digits=2)
@
This output is similar to the gene-level summaries reported by
\verb=skatMeta=, but also includes additional information.  Though we employ a score test, which does not explicitly estimate
genetic effects, we do report estimated effects \verb=beta=
and their standard error \verb=se=.  These can be thought of as one-step approximations to
standard maximum likelihood estimates, which
may differ slightly when effect sizes are very large
\citep{voorman2012fast}.  The additional columns \verb=cmafTotal= and
\verb=cmafUsed=, \verb=nsnpsTotal= and \verb=nsnpsUsed= distinguishing the total cumulative minor allele frequency and number of
SNPs, from those that are used in the test.  Also, note the genes for which no p-value is returned.  In these genes, no SNPs met the
inclusion criteria, and thus burden is zero for all individuals, and
the test is undefined.

\subsection{Meta analyzing SKAT-O}\label{sec:skato}
Recently, the `optimal' SKAT \citep[SKAT-O, ][]{lee2012optimal} was proposed to
test weighted averages of SKAT and burden tests of the form
\begin{align*}
Q_o(\rho) &= (1-\rho) \left(\sum_{j=1}^pw^{skat}_jU_j^2 \right)+ \rho
\left(\sum_{j=1}^p w^{burden}_jU_j\right)^2 \\
\end{align*}
When $\rho=0$ this gives the SKAT test $Q_o =\sum_{j=1}^pw_jU_j^2$ , and when
$\rho=1$ it gives the burden test $Q_o =
\left(\sum_{j=1}^p w_jU_j\right)^2$. Note that unlike the version of
SKAT-O implemented in the \pkg{SKAT} package, we do not assume that the weights used in the
burden test are the same as those used in SKAT.

When $\rho$ is fixed, computation of a p-value can be obtained
as with SKAT, using that the test statistic $Q_o(\rho)$ is a quadratic
form of normally distributed score vector $U$. \citet{lee2012optimal}  propose choosing the
value of $\rho$ which minimizes the $p$-value of the test over a set
values of $\rho$ in the interval $[0,1]$. That is, they find the
`optimal' linear combination of SKAT and burden tests. However, when $\rho$ is
chosen to minimize the $p$-value, the actual $p$-value reported must reflect
the flexibility afforded by this choice. For example, in the simple case where we
allow $\rho$ to be either $0$ or $1$, we would perform both a
burden test and a SKAT test and record the minimum of the
$p$-values. We might correct for multiple testing using a
 \v{S}id\'ak or Bonferroni correction.  However, this is not optimal, since burden tests and SKAT
are correlated with each other, especially when the number of variants
in a region is small.  As shown by \citet{lee2012optimal} it is possible to calculate the
distribution of the minimum $p$-value over any sequence of $\rho$'s,
accounting for this correlation between tests.
The mathematics are given in the Appendix of that paper, which require
individual level information.  In Section~\ref{sec:method} we
give an alternate derivation of the calculations which can be
performed at the meta-analysis level.

The syntax of \verb=skatOMeta= is similar to that of \verb=skatMeta=
and \verb=burdenMeta=, with the only difference being that weights for
SKAT and the burden test must be distinguished, and values of $\rho$
must be specified. These are given in the arguments

\begin{itemize}
 \item{\verb=skat.wts= and \verb=burden.wts=} Either a function, or a
   character string specifying a column in the SNPInfo file, which
   gives the weights to be used in SKAT and the burden test,
   respectively.  The default is to use the `beta' weights in SKAT and
   T1 weights for the burden test.
 \item{\verb=rho=} The values of $\rho$ to be used in SKAT-O. The
   default is \verb=c(0,1)=, which computes SKAT and the burden test,
   and reports the minimum $p$-value adjusted for multiple testing.
 \end{itemize}

 Here, we illustrate using SKAT-O using`Wu' weights in both the
 SKAT and the burden test, and choose a sequence of 11 $\rho$'s in $[0,1]$

<<>>=
meta.skato.results <- skatOMeta(c1, c2, rho=seq(0,1,length=11),
     burden.wts = function(maf){dbeta(maf,1,25)}, SNPInfo = SNPInfo, method = "int")

format(head(meta.skato.results),digits=2)
@ 

The format of the results is similar to the other region-based tests, with the
addition of the fields \verb=pmin=, which specifies the minimum $p$-value among
the tests considered, and \verb=rho=, which gives the value of $\rho$
which resulted in \verb=pmin=. Note that $p$ should always be at least as
large as \verb=pmin=, since SKAT-O corrects for the multiple tests
considered. If there is a single SNP in a gene, or either
SKAT or the burden test is undefined, \verb=pmin= will be identical to \verb=p=.

The additional column \verb=errflag= indicates a possibly inaccurate
$p$-values. A value of $0$ indicates no error, values
larger than $0$ indicate potentially inaccurate $p$-values. Genes where \verb=errflag= is
not zero can be re-run with more accurate calculation method, such as
\verb=saddlepoint=, which is slower, but more accurate,
than \verb=integration=. 

<<>>=
table(meta.skato.results$rho)
@ 

The table above displays the value of $\rho$ which gave the smallest
$p$-value among the 11 combinations of SKAT and the burden test. Note that the smallest $p$-value is
typically given by either SKAT or the burden test ($\rho=0$ or $\rho=1$), rather than a
proper combination of them.  For this reason, the default is simply
\verb=rho=c(0,1)= which performs only SKAT and the burden test, and
is typically much faster. 

Though \verb=skatOMeta= in some sense
supersedes \verb=skatMeta= and \verb=burdenMeta=, we
recommended running SKAT and burden tests on their own for
quality control. For instance, here we calculate SKAT and the burden
test that are used SKAT-O, and plot the SKAT-O
$p$-value, to the minimum $p$-value from SKAT and the burden test.

<<label=fig1plot,include=FALSE>>=
wu.burden <- burdenMeta(c1, c2, wts = function(maf){dbeta(maf,1,25)}, 
                        SNPInfo=SNPInfo)
pseq <- seq(0,1,length=100)
plot(y=meta.skato.results$p, x=pmin(wu.burden$p,meta.results$p), 
     xlab ="Minimum of SKAT and Burden", ylab = "SKAT-O")
abline(0,1,lty=2)
lines(x=pseq,y=1-(1-pseq)^2,col=2,lty=2,lwd=2)
legend("bottomright", lwd=2,lty=2,col=2,legend="Sidak correction")	
@


\begin{figure}
\begin{center}
<<label=fig1,fig=TRUE,echo=FALSE>>=
<<fig1plot>>
@
\end{center}
\caption{Comparison between SKAT-O $p$-values, and the minimum of SKAT
  and burden test $p$-values. The  \v{S}id\'ak correction for two
  independent tests is $1-(1-p)^2$.}
\label{fig:one}
\end{figure}

We see that the SKAT-O $p$-values are always larger than the
minimum of the SKAT and burden test $p$-values (the black line), but smaller than the
 \v{S}id\'ak correction for two independent tests (the red line). In addition to the \verb=errflag=
column, plots like these provide a rough check of the SKAT-O
$p$-value accuracy, which can occasionally be hard to achieve for
small $p$-values.

We also note here that the $p$-values for SKAT-O are close to \v{S}id\'ak correction for the minimum of SKAT and burden tests. As the number of variants increases, the SKAT and burden tests become independent. Thus, in the relatively large genes in this example, little efficiency is gained using the SKAT-O correction relative to the \v{S}id\'ak correction.

\subsection{Meta analyzing single SNPs}\label{sec:ss}

While \verb=skatMeta= and \verb=burdenMeta= perform region-based
tests, the function \verb=singlesnpMeta= can be used to perform score
tests for single SNP associations. The only additional option
available is whether or not to report cohort-specific effects and
standard errors (argument \verb=cohortBetas=)
<<>>=
meta.ss.results <- singlesnpMeta(c1, c2, SNPInfo = SNPInfo,
                                 cohortBetas = TRUE)
format(head(meta.ss.results),digits=2)
@ 
The output is similar to that of \verb=burdenMeta= and
\verb=skatMeta=, but now includes effect estimates for single cohorts,
suffixed by the name given to the corresponding \verb=skatCohort= object.

As with the burden test, the coefficients and standard errors can be thought of as
one-step approximations to maximum likelihood estimates.

\section{Binary and survival outcomes}\label{sec:other}
The \pkg{skatMeta} package can also handle binary data and survival
outcomes for unrelated
populations.  

For binary outcomes, this is specified with the \verb=family=
argument in the function \verb=skatCohort=, which accepts a family
object as in the function \verb=glm=.
<<>>=
c1.bin <- skatCohort(Z2, ybin~bmi+sex, family = binomial(), 
                     SNPInfo = SNPInfo, data = pheno1)
@

For survival outcomes this fitting process is somewhat more involved,
and handled with a separate function \verb=skatCoxCohort=.  The
\verb=formula= syntax is the same as that of \verb=coxph=:
<<>>=
c1.cox <- skatCoxCohort(Z=Z1, Surv(time, status) ~ bmi + strata(sex), 
                                 SNPInfo = SNPInfo, data =pheno1)
@ 
Regardless of model choice, the scores are meta-analyzed in exactly
the same way; no modification is necessary of the
\verb=skatMeta= function:
<<>>=
cohort1.bin.results <- skatMeta(c1.bin, SNPInfo = SNPInfo, 
                         aggregateBy = "gene")
cohort1.cox.results <- skatMeta(c1.cox, SNPInfo = SNPInfo, 
                         aggregateBy = "gene")
@ 

\section{Filtering}\label{sec:filter}
Sometimes it will be desirable to remove SNPs from the analysis, based
on minor allele cutoffs or other criteria. Some `static' criteria, such as SNP
annotation, can be handled by subsetting the SNP Info file given to
\verb=skatMeta=, \verb=burdenMeta= or \verb=singlesnpMeta=.

A common filtering criteria is minor
allele frequency, which in general depend on the individuals
contributing to an analysis. For this reason,\verb=skatMeta=,
\verb=burdenMeta= or \verb=singlesnpMeta= have an
option \verb=mafRange= which includes SNPs within a specified range of minor
allele frequencies.  For instance, if one only wishes to use only 
variants with MAF between 0 and 0.05, one would perform the previous
analysis by
<<>>=
meta.results <- skatMeta(c1, c2, SNPInfo = SNPInfo, 
                         aggregateBy = "gene", mafRange = c(0,0.05))
head(meta.results)
@
Note that the number of SNPs now included is lower than the previous analysis.

\section{Conditional analyses}\label{sec:conditional}
In many cases, it will be desirable to adjust for SNPs in a
conditional analysis. This can be done in the package with the functions
\verb=skatCohortAdj= and \verb=skatFamCohortAdj=. Note that these adjustments
are required by the individual cohorts, and are not performed at the
meta-analysis. The syntax and
output are identical to \verb=skatCohort= and \verb=skatFamCohort=,
but require the additional argument \verb=adjustments=. The
\verb=adjustments= should be in the same format as the SNPInfo file
(i.e. containing a `gene' and `Name' column specifying SNP names, or
designated alternative columns), and specifies which genes to adjust
for which SNPs. For instance, in the below example, we condition
`gene1' analyses on the SNPs labeled 1000001, 1000002, and 1000003,
`gene2' on SNP 1000020, and `gene13' on SNP 1000100.

<<>>=
adjustments <- SNPInfo[c(1:3, 20,100), ]
adjustments

####run on each cohort:
c1.adj <- skatCohortAdj(Z=Z1, y~sex+bmi, SNPInfo = SNPInfo, 
        adjustments=adjustments, data =pheno1)
c2.adj <- skatFamCohortAdj(Z=Z2, y~sex+bmi, SNPInfo = SNPInfo, 
        adjustments=adjustments, fullkins=kins, data=pheno2)

SNPInfo.sub <- subset(SNPInfo, (SNPInfo$gene %in% adjustments$gene) & 
        !(SNPInfo$Name %in% adjustments$Name) )

#skat
out.skat <- skatMeta(c1.adj,c2.adj, SNPInfo = SNPInfo.sub)
head(out.skat)
@ 

The output of \verb=skatCohortAdj= and \verb=skatFamCohortAdj= have class
\verb=skatCohort=, and contain only the genes listed in the
\verb=adjustments= argument. Any of the meta-analysis functions can be
applied to these objects as before.  In the above example, we
performed SKAT, using the subset of the full SNPInfo file relevant to
the conditional analyses.

\section{Parallelization}
\label{sec:bigdata}
In most cases, results for individual cohorts can be computed for
whole-exome and exome-chip on several thousand subjects in a few
minutes, without use of specialized data structures or parallel
processing. However, if either memory limitations or processing speed make this
cumbersome, we provide a generic concatenation function
\verb=c(...)= to combine
multiple \verb=skatCohort= objects. This can be employed to break up
computation over, e.g. chromosomes. For instance, we can separate
computation for the first cohort over the first and second set of 50
genes.  We perform this example in sequence, but it would be a simple
matter to run each subset of genes in parallel.

<<>>=
##subset SNPInfo file to first 50 genes, and second 50 genes:
SNPInfo1 <- subset(SNPInfo, SNPInfo$gene %in% unique(SNPInfo$gene)[1:50] )
SNPInfo2 <- subset(SNPInfo, !(SNPInfo$gene %in% unique(SNPInfo$gene)[1:50]) )

##subset corresponding genotype files:
Z1.1 <- subset(Z1, select = colnames(Z1) %in% SNPInfo1$Name) 
Z1.2 <- subset(Z1, select = colnames(Z1) %in% SNPInfo2$Name) 

##run skatCohort separately on each chunk:
c1.1 <- skatCohort(Z1.1, y~bmi+sex, SNPInfo = SNPInfo1, data = pheno1)
c1.2 <- skatCohort(Z1.2, y~bmi+sex, SNPInfo = SNPInfo2, data = pheno1)

##combine results:
c1 <- c(c1.1, c1.2)
class(c1)
@ 

We do not provide a method for computing the meta analysis in
parallel.  Individual \verb=skatCohort= objects are usually under 20
megabytes, and can be processed sufficiently quickly that
parallelization is typically not warranted.

\section{Method details}
\label{sec:method}
In this section we describe the method used in the meta-analysis. The theory for
meta-analyzing score tests can be found in e.g. \citet{lin2010relative}. 

We first focus on continuous and binary outcomes; modifications for
family data and survival outcomes are described subsequently.  For cohort $k$, denote $y^{(k)}$ as the $n$-vector of outcomes,
$\hat y^{(k)}$ as the corresponding vector of fitted values under the
null-model (i.e. without
genotype), and $G^{(k)}$ as the $n \times p$ matrix of genotypes in
one gene, and $X^{(k)}$ as the
$n \times q$ matrix of adjustment variables.

\subsection{Cohort level functions}
For each gene under
consideration, \verb=skatCohort= computes the vector of scores $U^{(k)}$
and their corresponding variances $V^{(k)}$.  For unrelated
individuals these are given by
\begin{align*}
  U^{(k)} & = G^{(k)T}(y^{(k)}-\hat y^{(k)})/\sigma^2_k \in \mathbb{R}^p \\
  V^{(k)} & = (W^{(k)} G^{(k)})^T (I-H^{(k)}) W^{(k)} G^{(k)}/\sigma^2_k \in
  \mathbb{R}^{p \times p},
\end{align*}
where $\sigma^2_k$ is the the residual variance for continuous data
and is 1 for binary data, 
$$
W^{(k)} = 
\begin{cases} 
\mathrm{diag}(\sqrt{\hat y^{(k)}(1-\hat y^{(k)})}) & \text{for
  binary data} \\
I_n & \text{for
  continuous data} \\
\end{cases},
$$
and 
$$H^{(k)} =  X^{(k)}W^{(k)} ( (W^{(k)}X^{(k)})^T W^{(k)}X^{(k)})^{-1}( W^{(k)}X^{(k)})^T,$$
is a projection matrix. These can easily be obtained from the output
of the \verb=glm= function.

The \verb=skatCohort= function computes these for each gene under
consideration, and concatenates them in a list.  If genotype is missing, it is imputed to the
cohort-specific minor allele frequency. This is done one gene at a
time as scores and information are computed, rather than modifying the
genotype argument \verb=Z=, which prevents \proglang{R}! from copying
potentially large genotype matrices.  In addition, the
sample size and minor allele count are also stored. 

\subsection{Meta analysis}
Each function \verb=skatMeta=, \verb=burdenMeta=, and \verb=singlesnpMeta=
first begins by meta-analyzing the scores and their variances across
all cohorts, for each unique gene name in the \verb=SNPInfo= file. This follows the simple formulas  
$$
U = \sum_k U^{(k)},\quad 
V =  \sum_k V^{(k)}.
$$
Under the null hypothesis of no genetic effect, we have that
$$V^{-1/2}U \rightarrow_d N_p(0, I_p)$$
or that $U$ is approximately $N_p(0,V)$

Each of the tests statistics computed computed (SKAT, burden, and
single snp) are functions of the vector $U$. 
\paragraph{skatMeta}
For a vector of weights $w=(w_1, \dots, w_p)^T$ denote $R= \mathrm{diag}(\sqrt{w_1}, \dots, \sqrt{w_p})$
The SKAT statistic $Q$ is given by
$$ Q = \sum_{j=1}^p w_j U_j^2 = \|RU\|^2_2,$$
where $S_j = \sum_k {g^{(k)}_j}^T(y^{(k)}-\hat y^{(k)})/\sigma^{2}_k$
and $g_j^{(k)}$, is the genotype for SNP $j$ in cohort $k$. Asymptotically
$$Q \sim\sum \lambda_i \chi^2_1,$$ where the $\lambda_i$ are the
eigenvalues of $RVR$.  The distribution function of a sum of eigenvalues can be
approximated using the function \verb=pchisqsum= in the \pkg{survey}
package \citep{lumley2004analysis}.  By default, we use \verb!method= 'saddlepoint'! for the
saddle point approximation, which is implemented in pure \proglang{R}.  A
slightly faster option is the Davies method, which inverts the
characteristic function, and requires the \pkg{CompQuadForm} package \citep{compquadform}.

\paragraph{burdenMeta and singlesnpMeta} The computation for burden tests is much more
straightforward. For a weight vector $w=(w_1, \dots, w_p)^T$, the
vector of burdens is given by $G^{(k)}w$, and the score for its association
by $w^TG^{(k)T}(y^{(k)}-\hat y^{(k)})/\sigma^2_k$. Summing the scores
across cohorts, we get the test statistic
$$T = \sum_k w^TG^{(k)T}(y^{(k)}-\hat y^{(k)})/\sigma^2_k = w^TU$$ 
which is approximately $N_1(0,w^TVw)$. The function
\verb=singlesnpMeta= is a special case of this, where $w_j=1$ for the
variant of interest, and is otherwise zero.

In order to approximate the coefficients and standard errors, we note
that the first step of the Fisher scoring algorithm estimates the
coefficients and their standard errors from the first and second derivatives of
the likelihood, which are given by the score and its variance.  We
plug in the meta-analyzed score and variance to obtain the one-step approximations
$$\hat \beta = \frac{w^TU}{w^TVw}, \text{ and }\quad \hat{se}(\hat\beta)= \frac{1}{\sqrt{w^TVw}}.$$

If we apply the Wald test to these coefficients, the resulting
$\chi^2$ test statistic is exactly the $\chi^2$ test statistic for the score test.


\paragraph{SKAT-O}
Denote the SKAT-O statistics
$$Q_o(\rho) = (1-\rho)\|UR\|_2^2 + \rho (w^TU).$$
Our goal is to find the value of $\rho$ which gives the smallest
$p$-value among a sequence of $\rho$'s (given in \verb=rho=), and correct for the
flexibility afforded by this choice. 

Denote the sequence of $\rho$'s $0 \leq \rho_1<
\dots<\rho_q\leq 1$ at which we calculate
$p$-values for $Q_o(\rho)$. In the simplest case, we can choose
$\rho_1=0$ and $\rho_2 =1$ to perform the SKAT and burden tests. We
then calculate the $p$-values for each of these tests. For fixed $\rho$ , $ Q_o(\rho_i)$ follows a mixture of $\chi^2_1$
variables, where the mixing parameters are given by the eigenvalues
of $W^{1/2}VW^{1/2}$, where $W = ((1-\rho) \cdot RR +
\rho \cdot ww^T)$. To calculate $W^{1/2}$ we use the eigenvalue
decomposition of the matrix  $((1-\rho) \cdot RR +
\rho \cdot ww^T)$. 

Now, let $p_i$ be the $p$-value for $Q_o(\rho_i)$. We then select the minimum $p$-value $p_{min} =
\min\{p_1,\dots,p_q\}$. To find the actual $p$-value we wish to
report, we must compute
the probability that of observing a $p_{min}$ smaller than what was
observed under the null hypothesis. Unfortunately, it is difficult to directly
write down a distribution of $p_{min}$ under the null hypothesis using
well-known distributions. However, the distribution of $p_{min}$ is
numerically tractable if we break the distribution of $Q_o(\cdot)$
into the conditional distribution of the SKAT statistic given the
value of the burden test, and then average over values of the burden
test.  In order to do this, first note that if $T_i$ is the $1-p_{min}$ quantifier of the distribution of
$Q_o(\rho_i)$, we know that the test using $\rho_i$ will give a
$p$-value less than $p_{min}$ when $(1-\rho_i)\|UR\|_2^2 +\rho_i
(w^TU)^2 > T_i$, i.e. when 
$$\|UR\|_2^2 > \frac{T_i - \rho_i (w^TU)^2}{1-\rho_i}.$$
Thus, given the value of the burden test, we observe a smaller
$p$-value than $p_{min}$ when
$$\|UR\|_2^2 > \min_{i=1\dots q} \left\{ \frac{T_i - \rho_i (w^TU)^2}{1-\rho_i}\right\}.$$
So, it suffices to find the conditional distribution of $UR$ given
$w^TU$, calculate $Pr(\|UR\|_2^2 > \min_{i=1\dots q} \left\{ \frac{T_i
    - \rho_i (w^TU)^2}{1-\rho_i}\right\})$, and integrate it over the
distribution of $w^TU$.  Since $UR$ and
$w^TU$ are linear combinations of a Normal vector, we can write down
their joint distribution.  Standard algebra gives that 
$$UR\mid\{ w^TU = a \} \sim N_p\left(\frac{RVw}{w^TVw}a, \quad R^TVR - \frac{RVww^TVR}{w^TVw}\right).$$
We can again determine the distribution of $\|UR\|_2^2$ using normal
quadratic forms. These differ slightly from the calculations used in
SKAT, due to the presence of the mean vector
$\frac{RVw}{w^TVw}a$. Instead, the distribution of $\|UR\|_2^2 \mid \{
w^TU = a \}$ is the sum of non-central $\chi^2$ distributions, whose
distribution functions are available in any of the functions in the
\pkg{CompQuadForm} package \citep{compquadform}. For speed, we recommend \verb!method=`integration'!,
which first attempts the Davies method, and in the case where this gives an error, it uses the
\verb=farebrother= method. Using \verb!method=`saddlepoint'! is slower, but has
higher relative accuracy.

Thus, if we denote $\phi(\cdot)$ as the marginal
distribution of $w^TU$, which is $N(0, w^TVw)$, we can find the true significance of $p_{min}$
using the following algorithm
\begin{enumerate}
 \item For each $\rho_i$, calculate $p_i$ based on the quadratic forms $Q_o(\rho_i)$. 
   \item Calculate $p_{min} = \min\{p_1,\dots,p_q\}$. 
     \item Calculate $T_i$'s, the $(1-p_{min})$ quantiles of the
       $Q_o(\rho_i)$.  Denoting $F_{Q_i}$ as the distribution function of
       $Q_o(\rho_i)$ found in Step 1, we do this by solving
       $F_{Q_i}(x) =p_{min}$ for $x$ using \verb=uniroot=. Potential inaccuracies here
     result in \verb!errflag=2!.
       \item Form the conditional distribution  $UR\mid \{ w^TU=
         a\}$, and the conditional probability $$f_S(a) = Pr\left(\|UR\|_2^2 > \min_{i=1\dots q} \left\{ \frac{T_i
    - \rho_i (w^TU)^2 }{1-\rho_i} \right\} \mid w^TU=
         a\right).$$
       \item Numerically integrate
         $$p_{actual} = \int_{-\infty}^{\infty} f_{S}(a)\phi (a)da$$
         and report $p_{actual}$.
   \end{enumerate}


\subsection{Family data}
This framework extends easily to the case of family data.  A more detailed description of this method is given in \citet{chen2012sequence}. In the pooled-data analysis, the related individuals would have a different form of the likelihood, corresponding to a mixed-model.  But, since likelihoods are on the same scale, the scores from related and un-related data can be combined as usual. We can think of all cohorts using a mixed-model, but where un-related individuals have no random effects. 

Here, instead of outcomes $y$ having variance $\sigma^2_y I_n$, they
have variance $\sigma^2_y\Omega =  \sigma_y^2(\theta_1\Theta +
\theta_2 I_n)$, where $\Theta$ is twice the kinship matrix, and $\theta_1 + \theta_2 =1$. In unrelated individuals $\Theta = I_n$.  The scores for this model are
$$U^{(k)} = {G^{(k)^T}} \Omega^{-1} (y^{(k)}-\hat y^{(k)}).$$
If the matrix of covariates is $X$, the variance of the score is
$$
V^{(k)} = G^{(k)T}\left(\Omega^{-1} - \Omega^{-1}X(X^T\Omega^{-1}X)^{-1}X^T\Omega^{-1}\right)G^{(k)}.
$$
We calculate these using the \verb=lmekin= function in
the \pkg{coxme} package \citep{coxme}. These can be combined as usual with scores and variances from other
cohorts. 

Unfortunately, the score test for generalized linear mixed models does
not exist in closed form, and no available implementation of
generalized linear mixed models allows large used-defined variance
components, such as a kinship matrix.  We therefore leave the
extension of binary traits for family data for future work.

\subsection{Survival outcomes}
We apply a similar procedure to obtain test statistics for Cox
proportional hazards regression. In this case, the score test is
anti-conservative, so we instead use a signed likelihood ratio
statistic \cite{lumley2013partial}. Specifically, let $S^2_j$ is the likelihood ratio test
statistic for the including of the single variant $j$ in the null
model, $\beta_j$ be the corresponding maximum likelihood estimate, and
$V$ be the variance of the usual score test. 
we use $$U_j = \text{sign}(\beta_j) S_j/\sqrt{V_{jj}},$$
which is asymptotically equivalent to the score test, but enjoys
superior finite sample performance. 

In order to implement this procedure efficiently, we use a
modified version of the\verb=coxph.fit= function in
the \pkg{survival} package \citep{survival}.



\section{Summary}
In this paper, we have described the basic usage and syntax of the
\pkg{skatMeta} \proglang{R} package.  This includes functions for computing
cohort level results with the functions \verb=skatCohort= for
unrelated individuals, \verb=skatFamCohort= for related individuals,
and \verb=skatCoxCohort= for survival outcomes. We also describe the
functions that pool cohort level results for meta analysis, including
\verb=skatMeta=, \verb=burdenMeta=, \verb=skatOMeta= and \verb=singlesnpMeta=.  
% where $W = \mathrm{diag}(w_1,\dots,w_p)$, $G^{(k)} = [g^{(k)}_1,\dots, g^{(k)}_p]$ and $P^{(k)}$ is the adjustment matrix for covariates (i.e. $P^{(k)}y^{(k)} = y^{(k)}-\hat y^{(k)}$).  

%The quantities which we can meta-analyze are the scores
%$$ \frac{1}{\sqrt{n}} \sum_k \frac{1}{\sigma^{2}_{y,k}} {g^{(k)}_j}^T(y^{(k)}-\hat y^{(k)}) \qquad j = 1\dots p,$$
%and the score variance, or the Information
%$$\frac{1}{n}\sum_k \frac{1}{\sigma^{2}_{y,k}} {G^{(k)}}^TP^{(k)}G^{(k)}.$$
%It is easy to see that the quantities ${g^{(k)}_j}^T(y^{(k)}-\hat y^{(k)})$ and ${G^{(k)}}^TP^{(k)}G^{(k)}$ can be calculated in each cohort separately. That is, we can computed the SKAT statistic \emph{exactly} as it would if all data were pooled together\footnote{assuming covariates are fit separately in each cohort}. This is a general feature of score tests, as noted by \citet{lin2010relative}.

\bibliographystyle{plainnat}

\input{skatMeta.bbl}

\end{document}
