Two-Stage Testing in a high dimensional setting Jonker et al

Marianne Jonker, Radboudumc Nijmegen, the Netherlands

# Two-stage testing for interactions in a high dimensional setting

Marianne A Jonker  Luc van Schijndel  Eric Cator  Research Institute for
Medical Innovation, Science department IQ Health, Section Biostatistics,
Radboud University Medical Center, Nijmegen, Netherlands  Faculty of Science,
Mathematics Department, Radboud University, Nijmegen, Netherlands
<EMAIL_ADDRESS>
###### Abstract

[Abstract] In a high dimensional regression setting in which the number of
variables ($p$) is much larger than the sample size ($n$), the number of
possible two-way interactions between the variables is immense. If the number
of variables is in the order of one million, which is usually the case in
e.g., genetics, the number of two-way interactions is of the order one million
squared. In the pursuit of detecting two-way interactions, testing all pairs
for interactions one-by-one is computational unfeasible and the multiple
testing correction will be severe. In this paper we describe a two-stage
testing procedure consisting of a screening and an evaluation stage. It is
proven that, under some assumptions, the tests-statistics in the two stages
are asymptotically independent. As a result, multiplicity correction in the
second stage is only needed for the number of statistical tests that are
actually performed in that stage. This increases the power of the testing
procedure. Also, since the testing procedure in the first stage is
computational simple, the computational burden is lowered. Simulations have
been performed for multiple settings and regression models (generalized linear
models and Cox PH model) to study the performance of the two-stage testing
procedure. The results show type I error control and an increase in power
compared to the procedure in which the pairs are tested one-by-one.

KEYWORDS:
Epistasis, FWER, GWAS, independent stages, multiple testing correction,
bladder cancer

††articletype: Research Article

##  1 Introduction

Let, for an individual, $Y$ be an outcome of interest and
${{X}}^{t}=(x_{1},\ldots,x_{p})$ a vector of covariates. Our aim is to detect
covariates which are associated with the outcome $Y$ via an interaction. We
have independent data of $n$ individuals from a population of interest, with
$n\ll p$ (the number of individuals is much smaller than the number of
covariates). This research question is important in, for instance, genetics
and genomics, where the number of covariates are usually in the order of
magnitude of several millions.

In genetics, the conventional approach to searching for pairwise interactions
(between genetic markers) is to fit regression models and test for two-way
interactions for all covariate pairs one-by-one. The number of pairs grows
rapidly with the number of covariates. If data of one million covariates is
available, the number of covariate pairs is of the order of one million
squared. This is not only computational a burden, but the multiple testing
correction will wipe out all effects that are present in the data. A wide
variety of possible strategies have been proposed to search for interactions
between genetic markers and between markers and environmental factors (see
e.g., Kooperberg et al (2008), Murcray et al (2009), Gauderman et al (2010),
Dai et al (2012), Hsu et al (2012), Lewinger et al (2013), Pecanka et al
(2017), Wong (Chapter 6, 2021), Johnsen et al (2021), Kawaguchi et al (2022),
Kawaguchi et al (2023)). An interesting approach that tackles both the
computational problem and the multiplicity issue, is a two-stage approach. The
first stage is seen as a screening stage to identify and select covariates
that are more likely to be involved in an interaction. Only the selected
covariates are passed on to the second stage. Subsequently, in the second
stage (the verification stage), only pairs of covariates that were selected in
the first stage are tested for interaction.

By selecting a minor part of the covariates in the first stage, the
computational burden of the whole testing procedure is reduced provided that
the screening in the first stage is computationally simple. Even more
important, with a two-stage procedure the multiplicity correction in the
second stage may be reduced, which could result in a considerable power gain.
If the statistical procedures in the two stages are statistically independent,
the family wise error rate (FWER) of the whole testing procedure is controlled
if multiplicity correction in the second stage is applied for the tests that
are actually performed in the second stage. Independence between stages can be
obtained in multiple ways. The easiest way is to split the data, one part for
the first stage and the other part for the second stage. Another way to obtain
independence is “independence by design”; the test-statistics in both stages
use all data, but are, nevertheless, statistically independent. A third option
is to obtain independence by regressing out all the information in the data
that was used in the first stage from the test-statistics in the second stage
(see Pecanka et al (2017)). An overview of multiple two-stage testing
procedures with independent stages for a binary outcome is given by Wong, 2021
(Chapter 6, by Pecanka et al).

Kooperberg et al (2008) proposed a two-stage procedure for generalized linear
models (GLMs). In the first stage, for every covariate a regression model with
only the marginal effect is fitted and association between the covariate and
the outcome variable is tested. If the corresponding p-value is below a pre-
specified threshold, the covariate is passed on to the second stage. In the
second stage interaction is tested for every pair of selected covariates. This
is done by fitting a multivariable regression model and testing for an
interaction effect. Dai et al (2012) adjusted the procedure in order to detect
interactions between genetic markers and environmental factors. In the first
stage they only made a selection of the genetic markers and passed on the
environmental factors to the second stage without testing. They proved FWER
control of their procedure. In this paper, we focus on testing for
interactions between genetic markers in GLMs and the Cox model. Although the
proposed testing procedure is similar to the one in Dai et al (2012), proving
control of the FWER is more complex, because in Dai et al (2012) the
environmental factors are not screened in the first stage. We prove FWER
control, by showing asymptotic independence between the testing procedures in
the two stages. We also prove independence between the two stages (and thus
FWER control) for the linear regression model with unknown variance of the
error term (this proof of FWER control seems to be missing in Dai et al
(2012)). For the semi-parametric Cox proportional hazards model for time-to-
event outcomes, the proof of independence between stages relies on asymptotic
properties of the partial likelihood.

The remainder of the paper is organized as follows. In the sections 2 and 3
the theory of the two-stage testing procedure is explained for GLMs and the
Cox PH model, respectively. To study control of the FWER and the power gain,
multiple simulation studies have been performed for the linear, Poisson and
Cox models. This, and an application to bladder cancer data, are described in
Section 4. In Section 5 assumptions and further research is discussed. The
paper ends with multiple appendices with proofs of theorems and lemmas.

##  2 Two-stage testing procedure

###  2.1 Two-stages

To study the association between an outcome of interest $Y$ with an
interaction between variables $x_{k}$ and $x_{\ell}$, the data are usually
modeled with a GLM:

$\displaystyle
h\big{(}\mathbb{E}(Y|(x_{k},x_{\ell}),{\beta},{\nu})\big{)}=\beta_{0}+\beta_{1}x_{k}+\beta_{2}x_{\ell}+\beta_{3}\;x_{k}x_{\ell},$
(1)

where $h$ is the (chosen) link function and $\beta_{0},\ldots,\beta_{3}$ are
the unknown regression parameters of which $\beta_{3}$ reflects the
interaction effect. The unknown parameter ${\nu}$ is an unknown nuisance
parameters, e.g., the variance of the error term in the linear regression
model. The parameters can be estimated by their maximum likelihood estimators
and the null hypothesis $H_{0}\mathrel{\mathop{\mathchar
58\relax}}\beta_{3}=0$ tested with, for instance, a Wald test.

In a high dimensional setting with many variables $x_{1},\ldots,x_{p}$
compared to the sample size $n$, $p\gg n$, it would not be sensible to test
all possible pairs of variables one-by-one, as the number of pairs grows
quadratically with the number of variables. Kooperberg et al (2008) proposed
the following two stage approach.

Stage 1, screening stage $S_{1}$
In the first stage all variables of ${{X}}^{t}=(x_{1},\ldots,x_{p})$ are
screened for association with the outcome $Y$ in a marginal model. More
specific, for variable $x_{k}$ the model

$\displaystyle
h\big{(}\mathbb{E}(Y|{{X}},{\gamma}=(\gamma_{0},\gamma_{1}),{\nu})\big{)}$
$\displaystyle=\gamma_{0}+\gamma_{1}x_{k}$ (2)

is fitted and the null hypothesis H${}_{0}\mathrel{\mathop{\mathchar
58\relax}}\gamma_{1}=0$ is tested with the Wald test-statistic
$T_{k}^{S_{1}}=\hat{\gamma}_{1}/(\widehat{\mbox{var}}\;\hat{\gamma}_{1})^{1/2}$
at a chosen significance level $\alpha_{1}$. Here $\hat{\gamma}_{1}$ is the
maximum likelihood estimator of $\gamma_{1}$ and
$\widehat{\mbox{var}}\;\hat{\gamma}_{1}$ an estimator of its corresponding
variance. If the null hypothesis is rejected, variable $x_{k}$ is passed on to
the second stage, the verification stage. This is done for every variable in
${X}$.

Stage 2, verification stage $S_{2}$
Interaction for all pairs of variables in stage 2 are tested. More specific,
suppose the variables $x_{k}$ and $x_{\ell}$ were passed on to the second
stage, then the model

$\displaystyle
h\big{(}\mathbb{E}(Y|{{X}},{\beta}=(\beta_{0},\ldots,\beta_{3}),{\nu})\big{)}=\beta_{0}+\beta_{1}x_{k}+\beta_{2}x_{\ell}+\beta_{3}x_{k}x_{\ell},$
(3)

is fitted and the null hypothesis H${}_{0}\mathrel{\mathop{\mathchar
58\relax}}\beta_{3}=0$ is tested with the Wald test-statistics
$T_{k\ell}^{S_{2}}=\hat{\beta}_{3}/(\widehat{\mbox{var}}\;\hat{\beta}_{3})^{1/2}$
at a chosen significance level $\alpha_{2}=\alpha/K_{1}$, with
$\hat{\beta}_{3}$ the maximum likelihood estimator of $\beta_{3}$, $K_{1}$ the
number of pairs in stage 2 and $\alpha$ the level of the test. This is done
for every pair of variables in stage 2.

With this two stage procedure, the computational burden decreases as the test
in the first stage is computational simpler than in the second stage (due to
dimension reduction of the parameter space) and interaction is tested for a
fraction of the variable pairs. Moreover, this approach is shown to be more
powerful than a single stage testing approach in which all pairs are tested
and a severe multiple testing correction is applied. This is discussed in the
next subsection.

We will consider a more general version of the two stage procedure; covariates
like gender, ethnicity or on behavior (e.g., smoking) that may be associated
with the outcome variable, but are not of interest for the interaction study
can be included in the regression models in both stages. Therefore, the models
in the two stages will be reformulated later.

###  2.2 Independent stages and multiple testing correction

To control the family wise error rate (FWER) of the two-stage testing
procedure, a multiple testing correction is applied in the second stage. This
multiple testing correction needs to take into account the fact that screening
has taken place in the first stage; viewing the data and screening for
interesting effects in the first stage cannot go unpunished. However, if under
the combined null hypothesis all test-statistics in the second stage are
statistically independent of those in the first stage, multiple testing
correction in the second stage is only needed for the number of tests that are
performed in that stage, as if no screening has taken place in stage 1. That
means that if $K_{1}$ pairs are tested in the verification stage, the
Bonferroni method will control the FWER of the whole procedure at a level
$\alpha$ by using $\alpha_{2}=\alpha/K_{1}$ as a corrected significance level
in the second stage. When the test-statistics in the two stages are not
statistically independent under the combined null hypothesis, the control of
the FWER is not guaranteed.111The number of pairs in the second stage,
$K_{1}$, is stochastic and a function of the test-statistics in the first
stage. This makes the two-stages dependent. We will come back to this
Subsection 2.4. Statistical independence between the two stages under the
combined null hypothesis is quite complex as “pairwise independence” is not
sufficient to control the type I error of the two stage testing procedure; we
need “complete independence”.

Pairwise independence. Suppose that for the variables $x_{k}$ and $x_{\ell}$
the model in (2) is fitted and the null hypotheses of no main effects are
tested with the test-statistics $T^{S_{1}}_{k}$ and $T^{S_{1}}_{\ell}$,
respectively. If both null hypotheses are rejected, it is tested whether the
interaction effect, $\beta_{3}$ in the model (3), equals zero with the test-
statistic $T^{S_{2}}_{k\ell}$. For pairwise independence the test-statistic
$T_{k\ell}^{S_{2}}$ must be independent of both $T_{k}^{S_{1}}$ and
$T_{\ell}^{S_{1}}$, for all $k$ and $\ell$.

Complete independence. Let ${\cal K}$ be the set of indices of the variables
that were passed on to the second stage. For complete independence, all test-
statistics in the second stage, $T_{k\ell}^{S_{2}},k,\ell\in{\cal K}$, should
be independent of all test-statistics $T_{j}^{S_{1}},j=1,\ldots,p$ in the
first stage jointly. Complete independence implies pairwise independence, but
the opposite is not true.

If pairwise and complete independence hold, the value of the significance
level in the first stage, $\alpha_{1}$, can be chosen freely. For high values
of $\alpha_{1}$ many variables will be passed on to the second stage and the
multiple testing correction in the second stage will be severe. The most
extreme case is $\alpha_{1}=1$ when no selection takes place in the first
stage. For a low value of $\alpha_{1}$ only a few variables will be passed on
and the multiple testing correction in the second stage is minor. Although
this may sound ideal, there exists the risk that variables that are associated
with the outcome via an interaction do not pass the first stage, which will
lower the power of the procedure. One may opt to use multiple thresholds in
the first stage. For variables which are a priori of more interest a higher
value for $\alpha_{1}$ (or even 1) can be taken than for variables that do not
have this preferred position.

###  2.3 Pairwise independence between stages

Before focusing on pairwise independence for the two stage testing procedure,
we formulate and prove pairwise independence in a more general setting. More
specific, the regression models are formulated without the interaction term,
but with non-specified extra covariates. These covariates can be chosen to be
equal to an interaction term. Consider the three nested generalized linear
regression models $M_{1},M_{2}$ and $M_{true}$:

$\displaystyle M_{1}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{X},{\beta}_{1},{\nu}_{1})\big{)}$
$\displaystyle=\beta_{1,1}x_{1}+\beta_{1,2}x_{2}+\ldots+\beta_{1,p_{1}}x_{p_{1}}$
(4) $\displaystyle M_{2}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{X},{\beta}_{2},{\nu}_{2})\big{)}$
$\displaystyle=\beta_{2,1}x_{1}+\beta_{2,2}x_{2}+\ldots+\beta_{2,p_{1}}x_{p_{1}}+\ldots+\beta_{2,p_{2}}x_{p_{2}}$
(5) $\displaystyle M_{true}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{X},{\beta}_{0},{\nu}_{0})\big{)}$
$\displaystyle=\beta_{0,1}x_{1}+\beta_{0,2}x_{2}+\ldots+\beta_{0,p_{1}}x_{p_{1}}+\ldots+\beta_{0,p_{2}}x_{p_{2}}+\ldots+\beta_{0,p_{0}}x_{p_{0}}$
(6)

with ${{X}}=(x_{1},x_{2},\ldots,x_{p})^{t}$, $h$ the link function, $p_{1}\leq
p_{2}\leq p_{0}\leq p$ and $x_{1}=1$ so that $\beta_{1,1},\beta_{2,1}$ and
$\beta_{0,1}$ are the intercepts in the models. The first index in
$\beta_{i,j}$ refers to the model and the second index to the covariate. Model
$M_{1}$ is nested in model $M_{2}$ and $M_{2}$ is nested in model $M_{true}$,
with the latter model assumed to be the true model. This can always be
effectuated, even if variables of the models $M_{1}$ and $M_{2}$ are not
associated with the outcome, by including variables in the true model (if they
are not in the model yet) with regression parameters equal to zero. Assume
that the conditional density $f_{{\theta}_{0}}$ of $Y|{X}$, distributed
according the true model $M_{true}$, belongs to a GLM with a canonical link
function and, therefore, can be written as:

$\displaystyle
f_{{\theta}_{0}}(Y|{X})=g(Y)\exp\big(V(Y)^{t}{Z}{\theta}_{0}-A({Z}{\theta}_{0})\big{missing}),$
(7)

where ${Z}$ is a function of the covariates ${X}$ and ${\theta}_{0}$ the
vector of unknown parameters in the true model; the regression and nuisance
parameters. The vector ${\theta}_{0}\in\mathbb{R}^{q_{0}}$ with $q_{0}-p_{0}$
the number of nuisance parameters in the model.

Let ${\theta}_{1}$ and ${\theta}_{2}$ be the parameter vectors for the models
$M_{1}$ and $M_{2}$ (i.e., regression and nuisance parameters). Their
parameter spaces are assumed to be subsets of $\mathbb{R}^{q_{1}}$ and
$\mathbb{R}^{q_{2}}$, respectively, with $q_{1}\geq p_{1}$ and $q_{2}\geq
p_{2}$ depending on the number of nuisance parameters. If no nuisance
parameters are present $\theta_{1},\theta_{2}$ and $\theta_{0}$ correspond to
the vector of regression parameters. The parameter spaces for the models
$M_{1},M_{2}$ and $M_{true}$ are nested. To project a parameter vector from
the space for $M_{1}$ or $M_{2}$ to the parameter space for model $M_{true}$
(the largest space), we define the projection matrices

$\displaystyle E_{1}$
$\displaystyle=\left(e_{1}|\ldots|\;e_{q_{1}}\right)\in\mathbb{R}^{q_{0}\times
q_{1}}$ $\displaystyle E_{2}$
$\displaystyle=\left(e_{1}|\ldots|\;e_{q_{2}}\right)\in\mathbb{R}^{q_{0}\times
q_{2}},$

in such a way that $\\{e_{1},\ldots,e_{q_{1}}\\}$ and
$\\{e_{1},\ldots,e_{q_{2}}\\}$ form orthonormal bases for the parameter spaces
for $M_{1}$ and $M_{2}$. To be more specific, these matrices represent linear
functions from the parameter spaces for $M_{i},i=1,2$ to $M_{true}$, which
will be used to project the regression coefficients within the two submodels
to the parameter space of the true model by left-multiplication. So, for
arbitrary parameter vectors ${\theta}_{1}\in\mathbb{R}^{q_{1}}$ and
${\theta}_{2}\in\mathbb{R}^{q_{2}}$, it holds that
$E_{1}{\theta}_{1}\in\mathbb{R}^{q_{0}}$ and
$E_{2}{\theta}_{2}\in\mathbb{R}^{q_{0}}$.

Let
${{X}}_{1}=(x_{1,1},\ldots,x_{1,p})^{t},\ldots,{{X}}_{n}=(x_{n,1},\ldots,x_{n,p})^{t}$
be a sample of variables and let $Y_{1}|{{X}}_{1},\ldots,Y_{n}|{{X}}_{n}$ be a
sample from the GLM $M_{true}$. The maximum likelihood estimators for
${\theta}_{1}$ and ${\theta}_{2}$ in the models $M_{1}$ and $M_{2}$ are
denoted as $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$. Furthermore, let
${\theta}_{1}^{0}$ and ${\theta}_{2}^{0}$ be the values that maximize the
pointwise limit of the log likelihood function corresponding to the different
models (divided by $n$). The projections of the vectors $\hat{\theta}_{2}$ and
${\theta}_{2}^{0}$ in the parameter space of model $M_{2}$ to the space that
is orthogonal to the parameter space of model $M_{1}$ are computed as
$\hat{\theta}_{2}^{\perp}=(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2})\hat{\theta}_{2}$
and
${\theta}_{2}^{0,\perp}=(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2}){\theta}_{2}^{0}$,
for $Id_{q}$ the $q\times q$ identity matrix. In the theorem below we will
prove that $\sqrt{n}(\hat{\theta}_{1}-{\theta}_{1}^{0})$ and
$\sqrt{n}(\hat{\theta}_{2}^{\perp}-{\theta}_{2}^{0,\perp})$ are asymptotically
independent if $M_{2}=M_{true}$. The proof of the theorem is given in Appendix
A.

######  Theorem 2.1.

Suppose that $Y_{1}|{{X}}_{1},\ldots,Y_{n}|{{X}}_{n}$ is a sample from the
generalized linear model with a canonical link function $M_{true}$ in (6). Let
$\hat{{\theta}}_{1}$ and $\hat{{\theta}}_{2}$ denote the maximum likelihood
estimators of ${\theta}_{1}^{0}$ and ${{\theta}}_{2}^{0}$ in the nested models
$M_{1}$ and $M_{2}$, respectively. If model $M_{2}$ coincides with the true
model $M_{true}$ (i.e., $p_{2}=p_{0}$), then
$\sqrt{n}(\hat{\theta}_{1}-{\theta}_{1}^{0})$ and
$\sqrt{n}(\hat{\theta}_{2}^{\perp}-{\theta}_{2}^{0,\perp})$ are asymptotically
independent.

Now, we translate the theorem to the two stage testing framework by
reformulating the GLMs so that it includes an interaction term and explicitly
defining the matrices $E_{1}$ and $E_{2}$. Consider the two nested GLMs
$M_{1}$ and $M_{2}$:

$\displaystyle M_{1}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{{X}},{\beta}_{1},{\nu}_{1})\big{)}$
$\displaystyle=\beta_{1,1}x_{1}+\beta_{1,2}x_{2}+\beta_{1,k}x_{k}$ (8)
$\displaystyle M_{2}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{{X}},{\beta}_{2},{\nu}_{2})\big{)}$
$\displaystyle=\beta_{2,1}x_{1}+\beta_{2,2}x_{2}+\beta_{2,k}x_{k}+\beta_{2,\ell}x_{\ell}+\beta_{2,k\ell}x_{k}x_{\ell}$
(9)

with ${{X}}=(x_{1},x_{2},\ldots,x_{p})^{t}$, $h$ the link function, $x_{1}=1$
so that $\beta_{1,1}$ and $\beta_{2,1}$ are the intercepts in the models, and
$x_{2}$ the “fixed” covariate which is included in all models and is not
tested (multiple “fixed” covariates can be included). The model $M_{1}$ equals
the GLM that is fitted in the first stage, the screening stage $S_{1}$, of the
two-stage testing procedure as defined earlier. In this stage it is studied
whether the covariate $x_{k}$ (genetic marker $k$) is associated with the
outcome $Y$ by testing the null hypothesis $H_{0}\mathrel{\mathop{\mathchar
58\relax}}\beta_{1,k}=0$. Model $M_{2}$ is the interaction model that is
fitted in stage 2, the verification stage $S_{2}$. Now, it is tested whether
the interaction term $x_{k}x_{\ell}$ is associated with $Y$, by testing the
null hypothesis $H_{0}\mathrel{\mathop{\mathchar 58\relax}}\beta_{2,k\ell}=0$.

Further, ${\beta}_{1}=(\beta_{1,1},\beta_{1,2},\beta_{1,k})^{t}$ and
${\beta}_{2}=(\beta_{2,1},\beta_{2,2},\beta_{2,k},\beta_{2,\ell},\beta_{2,k\ell})^{t}$
denote the vectors of regression parameters in the two models, and their
maximum likelihood estimators are denoted as $\hat{{\beta}}_{1}$ and
$\hat{{\beta}}_{2}$. The vectors
${\beta}_{1}^{0}=(\beta_{1,1}^{0},\beta_{1,2}^{0},\beta_{1,k}^{0})^{t}$ and
${\beta}_{2}^{0}=(\beta_{2,1}^{0},\beta_{2,2}^{0},\beta_{2,k}^{0},\beta_{2,\ell}^{0},\beta_{2,k\ell}^{0})^{t}$
are the values that maximize the pointwise limit of the log likelihood
function corresponding to the models $M_{1}$ and $M_{2}$ (divided by $n$).
These vectors coincide with the true parameters only if the corresponding
regression model equals the true model. Define
$\hat{\beta}_{2}^{\perp}=(\hat{\beta}_{2,\ell},\hat{\beta}_{2,k\ell})^{t}$ and
${\beta}_{2}^{0,\perp}=(\beta_{2,\ell}^{0},\beta_{2,k\ell}^{0})^{t}$, the
regression parameters for the covariates that are additional to model $M_{2}$
compared to model $M_{1}$.

If, in the model, there are no nuisance parameters (e.g., the Poisson and
logistic regression model), then ${\theta}_{1}={\beta}_{1}$ and
${\theta}_{2}={\beta}_{2}$. If, moreover, $\\{e_{1},\ldots,e_{q_{1}}\\}$ and
$\\{e_{1},\ldots,e_{q_{2}}\\}$ form canonical bases,
${\theta}_{2}^{0,\perp}=(0.\ldots,0,\theta_{q_{1}+1}^{0},\ldots,\theta_{q_{2}}^{0})^{t}$
and it follows from Theorem 2.1 that
$\sqrt{n}(\hat{\theta}_{1}-{\theta}_{1}^{0})=\sqrt{n}(\hat{\beta}_{1}-{\beta}_{1}^{0})$
and
$\sqrt{n}(\hat{\theta}_{2}^{\perp}-{\theta}_{2}^{0,\perp})=\sqrt{n}(\hat{\beta}_{2}^{\perp}-{\beta}_{2}^{0,\perp})$
are asymptotically independent. Thus,
$\sqrt{n}(\hat{\beta}_{1,k}-\beta^{0}_{1,k})$ and
$\sqrt{n}(\hat{\beta}_{2,k\ell}-\beta_{2,k\ell}^{0})$ are asymptotically
independent: the interaction effect $\hat{\beta}_{2,k\ell}$ in model $M_{2}$
(fitted in the second stage of the two-stage testing procedure) is
asymptotically independent of the marginal effects $\hat{\beta}_{1,k}$ in
model $M_{1}$ that is fitted in the first stage.

[Pairwise independence of estimated regression coefficients] Suppose that
$Y_{1}|{{X}}_{1},\ldots,Y_{n}|{{X}}_{n}$ is a sample from the GLM $M_{2}$ in
(9) with a canonical link function $h$. If there are no nuisance parameter,
$\sqrt{n}(\hat{\beta}_{1}-\beta_{1}^{0})$ and
$\sqrt{n}(\hat{\beta}_{2}^{\perp}-{\beta}_{2}^{0,\perp})$ are asymptotically
independent.

Suppose that in the GLMs with a canonical link function, the parameter vector
${\theta}_{0}$ as defined in the density of $Y|X$ in (7) is of the form
$(\nu_{0}^{t},{\beta}_{0}^{t})^{t}$, with $\nu_{0}$ a nuisance parameter. By
choosing $\\{e_{1},\ldots,e_{q_{1}}\\}$ and $\\{e_{1},\ldots,e_{q_{2}}\\}$ so
that they form canonical bases for the parameter spaces for $M_{1}$ and
$M_{2}$, Theorem 2.1 can be directly applied to prove pairwise independence;
i.e., asymptotically independence between the interaction effect
$\hat{\beta}_{2,k\ell}$ in the second stage and the marginal effects
$\hat{\beta}_{1,k}$ and $\hat{\beta}_{1,\ell}$ in the first stage of the two-
stage testing procedure. If there are nuisance parameters in the model, but
the parameter vector ${\theta}_{0}$ is not of the form
$(\nu_{0}^{t},{\beta}_{0}^{t})^{t}$ (the regression and the nuisance
parameters are not separated), like in the linear regression model (see
below), pairwise independence of the estimated regression coefficients is not
guaranteed by the theorem. As an example, suppose $Y|{X}$ has a normal
distribution with mean ${X}^{t}{\beta}_{0}$ and variance $\sigma_{0}^{2}$. By
writing the conditional density in the form (7):

$\displaystyle
f_{{\beta}_{0},\sigma_{0}^{2}}(y|{X})\;=\;\frac{1}{\sqrt{2\pi}}\exp\Big(y\;{X}^{t}{\beta}_{0}/\sigma_{0}^{2}-y^{2}/(2\sigma_{0}^{2})-({X}^{t}{\beta}_{0})^{2}/(2\sigma_{0}^{2})-\log\sigma_{0}\Big{missing})\;=\;\frac{1}{\sqrt{2\pi}}\exp\Big(V(y)^{t}Z\theta_{0}-A(Z\theta_{0})\Big{missing})$

it follow that the model belongs to the family of GLMs with a canonical link
function, with $V(y)=(y^{2},y)^{t}$ and
$A(Z{\theta}_{0})=({X}^{t}{\beta}_{0})^{2}/(2\sigma_{0}^{2})+\log\sigma_{0}$,
the parameter
${{\theta}_{0}}=(-1/(2\sigma_{0}^{2}),{{\beta}_{0}}^{t}/\sigma_{0}^{2})^{t}\in\mathbb{R}^{p+1}$,
and

$\displaystyle Z=\begin{pmatrix}~{}1&~{}~{}{\bf 0}^{t}~{}\\\
~{}0&~{}~{}{{X}}^{t}~{}\end{pmatrix}\in\mathbb{R}^{2}\times\mathbb{R}^{p+1}$
(10)

for ${\bf 0}\in\mathbb{R}^{p}$ a vector with zeroes. By Theorem 2.1, and for
the standard basis,
$\sqrt{n}(\hat{\beta}_{1,k}/\hat{\sigma}_{1}^{2}-\beta_{1,k}^{0}/\sigma_{1}^{2})$
and
$\sqrt{n}(\hat{\beta}_{2,k\ell}/\hat{\sigma}_{2}^{2}-\beta_{2,k\ell}^{0}/\sigma_{2}^{2})$
are asymptotically independent if model $M_{2}$ is the true model. Here
$\hat{\sigma}_{1}^{2}$ and $\hat{\sigma}_{2}^{2}$ are the estimated variances
of the errors in the models $M_{1}$ and $M_{2}$, $\sigma_{1}^{2}$ and
$\sigma_{2}^{2}$. From this independence results we can not (directly)
conclude asymptotic independence of
$\sqrt{n}(\hat{\beta}_{1,k}-\beta_{1,k}^{0})$ and
$\sqrt{n}(\hat{\beta}_{2,k\ell}-\beta_{2,k\ell}^{0})$. For the linear
regression model pairwise independence can be proven by straightforward
calculations.

######  Theorem 2.2 (Pairwise independence in the linear regression model).

Suppose that $Y_{1}|{{X}}_{1},\ldots,Y_{n}|{{X}}_{n}$ is a sample from the
linear regression model $M_{2}$ in (9) with the identical link function and
unknown variance $\sigma^{2}$. Then $\hat{\beta}_{1}$ and
$\hat{\beta}_{2}^{\perp}$ are independent.

The proof of Theorem 2.2 is given in Appendix A. In terms of the two-stage
testing procedure, the previous results state that if $M_{2}$ is the true
model, the estimate of the regression parameter for the interaction term in
the second stage, $\hat{\beta}_{2,k\ell}$ in the model $M_{2}$ in (9), is
(asymptotically) independent of the maximum likelihood estimators of the
regression parameter $\hat{\beta}_{1,k}$ in the model in (8) and by symmetry
also of $\hat{\beta}_{1,\ell}$.

The assumption that model $M_{2}$ is the true model is essential for the
proof. This assumption implies that all variables that are associated with the
outcome via a main effect or an interaction effect must be included in this
model. This assumption would be difficult to hold in a study in which one aims
to find new variables (interactions of genetic markers) that are associated
with the outcome. However, to control the FWER in the two-stage testing
procedure, independence is only required under the (full) null hypothesis.
Therefore, in the two-stage procedure just described, the null hypothesis
should state that there are no main and no interaction effects among the
variables that participate in the two-stage testing procedure to be associated
with the outcome; known main effects should be included as “fixed” covariates.
The paper by Dai et al (2012) may make you think differently about the choice
of the null hypothesis. However, in their paper pairwise independence between
the test-statistics in the two stages is proven under the (implicit)
assumption of no model misspecification of the interaction model in (9), and
thus under the assumption that all true main and interaction effects are
included in this model for every statistical test in the verification stage.
In other words, they (implicitly) work with the same null hypothesis. The
above discussion also means that, once we can prove complete independence, we
have FWER control in weak sense (and not in strong sense like is claimed by
Dai et al (2012)); under the union of all null hypotheses. A proof of Theorem
2.2 for unknown variance $\sigma^{2}$ seems to be missing in Dai et al (2012).

The theorems above state that there is pairwise independence between the
maximum likelihood estimators of the regression parameters. However, in the
two-stage testing procedure the test-statistic
$T_{k}^{S_{1}}=\hat{\beta}_{1,k}/\widehat{\mbox{var}}(\hat{\beta}_{1,k})^{1/2}$
for main effects in stage 1 has to be (asymptotically) independent of the
test-statistic
$T_{k\ell}^{S_{2}}=\hat{\beta}_{2,k\ell}/\widehat{\mbox{var}}(\hat{\beta}_{2,k\ell})^{1/2}$
for testing for interaction in stage 2, under the null hypothesis that all
main and interaction effects equal zero (so also
$\beta_{1,k}^{0}=\beta_{2,k}^{0}=\beta_{2,\ell}^{0}=\beta_{2,k\ell}^{0}=0$.
The Wald test-statistics in the first and second stage equal
$\hat{\beta}_{1,k}/\widehat{\mbox{var}}(\hat{\beta}_{1,k})^{1/2}=\sqrt{n}\hat{\beta}_{1,k}(I(\hat{{\theta}}_{1})^{-1}_{(k)})^{1/2}$
and
$\hat{\beta}_{2,k\ell}/\widehat{\mbox{var}}(\hat{\beta}_{2,k\ell})^{1/2}=\sqrt{n}\hat{\beta}_{2,k\ell}(I(\hat{{\theta}}_{2})^{-1}_{(k\ell)})^{1/2}$,
with $I(\hat{{\theta}}_{1})$ and $I(\hat{{\theta}}_{2})$ the estimated Fisher
information matrices in the models $M_{1}$ and $M_{2}$, and
$I(\hat{{\theta}}_{2})_{(k\ell)}^{-1}$ equal to the diagonal element of the
inverse Fisher Information matrix that corresponds to the parameter
$\beta_{k\ell}$ (and similar for $I(\hat{{\theta}}_{1})_{(k)}^{-1}$).
Asymptotic pairwise independence between the test-statistics follows from the
pairwise independence of the maximum likelihood estimators of the regression
parameters and Slutky’s lemma (see e.g., van der Vaart (1998)).

[Pairwise independence of test-statistics] Suppose that
$Y_{1}|{{X}}_{1},\ldots,Y_{n}|{{X}}_{n}$ is a sample from model $M_{2}$ in (9)
with a canonical link function $h$ without any nuisance parameter $\nu_{2}$,
or from the linear regression model (identical link function) with unknown
variance $\sigma^{2}$. Then, under the full null hypothesis of no main and
interaction effects, it holds that the Wald test-statistics $T_{k}^{S_{1}}$
and $T_{k\ell}^{S_{2}}$ are asymptotically independent.

###  2.4 Complete independence between stages

Define the generalized linear regression models:

$\displaystyle M_{1}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{X},{\beta}_{1},{\nu}_{1})\big{)}$
$\displaystyle=\beta_{1,1}x_{1}+\beta_{1,2}x_{2}+\beta_{1,j}x_{j}$ (11)
$\displaystyle M_{2}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{X},{\beta}_{2},{\nu}_{2})\big{)}$
$\displaystyle=\beta_{2,1}x_{1}+\beta_{2,2}x_{2}+\hskip
36.98866pt\beta_{2,k}x_{k}+\beta_{2,\ell}x_{\ell}+\beta_{2,k\ell}x_{k}x_{\ell}$
(12) $\displaystyle M_{true}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{X},{\beta}_{0},{\nu}_{0})\big{)}$
$\displaystyle=\beta_{0,1}x_{1}+\beta_{0,2}x_{2}$ (13)

where, this time, model $M_{1}$ (fitted in stage 1) is not nested in model
$M_{2}$ (fitted in stage 2), and the model $M_{true}$, the true model under
the full null hypothesis, is nested in both $M_{1}$ and $M_{2}$. Every model
has an intercept ($x_{1}=1$) and the “fixed” covariate $x_{2}$ (multiple
“fixed” covariates can be included).

Complete independence of the test-statistics can be proven under the null
hypothesis of no main and no interaction effects for independent
(standardized) covariates for the logistic, the Poisson and the linear
regression model. This will be proven in steps. The first step is mixed
pairwise independence. The proof of the theorem below is given in Appendix B.

######  Theorem 2.3 (Mixed pairwise independence of estimated regression
coefficients).

Suppose that $Y_{1}|{{X}}_{1},\ldots,Y_{n}|{{X}}_{n}$ is a sample from model
$M_{true}$ in (13) and let $\hat{{\beta}}_{1}$ and $\hat{{\beta}}_{2}$ be the
maximum likelihood estimators of the parameters
${\beta}_{1}=(\beta_{1,1},\beta_{1,2},\beta_{1,j})^{t}$ and
${\beta}_{2}=(\beta_{2,1},\beta_{2,2},\beta_{2,k},\beta_{2,\ell},\beta_{2,k\ell})^{t}$
in the models $M_{1}$ and $M_{2}$ in (11) and (12), respectively. Under the
null hypothesis of no main and no interaction effects (Model $M_{true}$) and
if the covariates $x_{j},x_{k},x_{\ell}$ are mutually independent, are
independent of the covariate $x_{2}$ and have mean zero, then
$\sqrt{n}(\hat{\beta}_{1}-{\beta}_{1}^{0})$ and
$\sqrt{n}\hat{\beta}_{2}^{\perp}$ are asymptotically independent in the
logistic, the Poisson and the linear regression model, with
$\hat{\beta}_{2}^{\perp}=(\hat{\beta}_{2,k},\hat{\beta}_{2,\ell},\hat{\beta}_{2,k\ell})^{t}$,
for all $j,k,\ell$ with $j\neq k,\ell$ and
$\hat{\beta}_{2}^{\perp}=(\hat{\beta}_{2,\ell},\hat{\beta}_{2,k\ell})^{t}$ if
$j=k$ (and similarly in case $k=\ell$).

In the theorem, ${\beta}_{1}^{0}=(\beta_{1,1}^{0},\beta_{1,2}^{0},0)$ with
$\beta_{1,1}^{0}=\beta_{0,1}^{0}$ and $\beta_{1,2}^{0}=\beta_{0,2}^{0}$, the
true parameter values in model $M_{true}$. From the theorem it follows that
$\sqrt{n}\hat{\beta}_{2,k\ell}$ and $\sqrt{n}\hat{\beta}_{1,j}$ are
asymptotically independent. If $j=k$ or $j=\ell$ the theorem implies pairwise
independence.

It directly follows from Theorem 2.3 and Slutky’s Lemma that the Wald test-
statistics in the second stage are asymptotically independent of the Wald
test-statistics in the first stage, under the null hypothesis of no main and
interaction effects (and under the assumptions stated in theorem above). This
leads to the following result.

[Mixed pairwise independence of test-statistics] Suppose that
$Y_{1}|{{X}}_{1},\ldots,Y_{n}|{{X}}_{n}$ is a sample from model $M_{true}$ in
(13). Under the null hypothesis of no main and no interaction effects (Model
$M_{true}$) and if the covariates $x_{j},x_{k},x_{\ell}$ are mutually
independent, independent of the covariate $x_{2}$ and have mean zero, then the
Wald test-statistics $T_{j}^{S_{1}}$ and $T_{k\ell}^{S_{2}}$ are
asymptotically independent in the logistic, the Poisson and the linear
regression model, for all $j,k,\ell$.

Complete independence is stronger than mixed pairwise independence of the
test-statistics. For complete independence every test-statistic
$T_{k\ell}^{S_{2}}$ in the second stage should be asymptotically independent
of all first stage test-statistics $T_{j}^{S_{1}},j=1,\ldots,p$ jointly.
Since, under the null hypothesis, the test-statistics in both stages are
asymptotically Gaussian, complete independence follows from the asymptotic
(mixed) pairwise independence.

Now, complete independence has been proven for the linear, logistic and the
Poisson regression models. For other GLMs complete independence may be true as
well, but this depends on the form of the model and needs to be checked for
each model separately. More details on how to check complete independence can
be found in Appendix B.

The number of pairs in the verification stage, $K_{1}$, determines the
multiple testing correction in stage 2. This number is stochastic and a
function of the test-statistics in the screening-stage. In order to obtain
independence between the two-stages, the number $K_{1}$ could be replaced by
the non-random value $(p\alpha_{1})^{2}$, which approximately equals the
expected number of pairs in the second stage under the (combined) null
hypotheses, so under a sparse setting.

Complete independence holds under the assumption that the covariates are
independent. In a GWAS the covariates (i.e., genetic markers) can not be
assumed to be independent; markers are expected to be correlated with
neighboring markers. This is only a fraction of the total number of markers in
the data set. In the simulation studies in Subsection 4.1 it will be seen that
in all considered settings with correlated markers the FWER is controlled.

Pairwise and complete independence imply FWER control of the proposed testing
procedure. However, under exactly the same assumptions as made in Theorem 2.3,
FWER control can also be proved directly and without replacing the stochastic
number $K_{1}$ in the multiple testing correction. This proof is given in
Appendix C. Summarized, the two-stage procedure in which a Bonferroni
correction is performed in the second stage only, preserves FWER control at
level $\alpha$ under the assumptions stated above.

##  3 Two-stage testing with the Cox proportional hazards model

The most popular model for time-to-event data is the Cox proportional hazards
(PH) regression model. This model is a semi-parametric model and does not
belong to the family of GLMs. The regression parameters in the model are
estimated by maximizing the partial likelihood. Asymptotic pairwise and
complete independence can be proved along the same line as before, but the
proof is based on the partial likelihood function (in stead of the full
likelihood).

Let $\tilde{T}$ be the time to an event of interest, from a well defined time
zero. Given a covariate vector ${X}$, the hazard function is defined as

$\displaystyle h(t|{X})=h_{0}(t)\exp({\beta}^{t}{X}),$

where $h_{0}$ is the completely unknown baseline hazard function and ${\beta}$
the vector of unknown regression coefficients. Let $C$ be the right-censoring
time, which is assumed to be independent of time $\tilde{T}$. Furthermore, let
$T=\min\\{\tilde{T},C\\}$ and $\Delta=1_{\\{\tilde{T}\leq C\\}}$ the indicator
function which indicates whether the event has been observed or is right-
censored. Suppose data of $n$ independent individuals is available:
$\\{(T_{1},\Delta_{1},{X}_{1}),\ldots,(T_{n},\Delta_{n},{X}_{n})\\}$, where
the index refers to the individual. In the Cox regression model the vector
with regression parameters ${\beta}$ is estimated by maximizing the partial
likelihood:

$\displaystyle
PL({\beta})=\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{\infty}\log\Bigg(\frac{\exp({\beta}^{t}{X}_{i})}{\sum_{j=1}^{n}Y_{j}(t)\exp({\beta}^{t}{X}_{j})}\Bigg{missing}){\rm
d}N_{i}(t),$ (14)

with $Y_{i}(t)=1_{\\{T_{i}>t\\}}$ indicating whether a person is still at risk
at time $t$ and $N_{i}(t)=1_{\\{T_{i}\leq t\\}}$ (see e.g., Aalen et al (2008)
and Therneau et al (2000)).

Let $M_{true}$ be the true model under the full null hypothesis, and let
$M_{1}$ and $M_{2}$ be the nested models with hazard functions

$\displaystyle M_{1}\mathrel{\mathop{\mathchar 58\relax}}\qquad h_{1}(t|{X})$
$\displaystyle=h_{0,1}(t)\exp(\beta_{1,1}x_{1}+\beta_{1,k}x_{k})$ (15)
$\displaystyle M_{2}\mathrel{\mathop{\mathchar 58\relax}}\qquad h_{2}(t|{X})$
$\displaystyle=h_{0,2}(t)\exp(\beta_{2,1}x_{1}+\beta_{2,k}x_{k}+\beta_{2,\ell}x_{\ell}+\beta_{2,k\ell}x_{k}x_{\ell}),$
(16) $\displaystyle M_{true}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h_{0}(t|{X})$ $\displaystyle=h_{0,0}(t)\exp(\beta_{0,1}x_{1}).$ (17)

Here, $x_{1}$ is the “fixed" covariate as explained before. The intercept is
missing, as it is absorbed in the baseline hazard function. The maximum
likelihood estimators of the regression parameters in the models $M_{1}$ and
$M_{2}$ are denoted as
$\hat{\beta}_{1}=(\hat{\beta}_{1,1},\hat{\beta}_{1,k})^{t}$ and
$\hat{\beta}_{2}=(\hat{\beta}_{2,1},\hat{\beta}_{2,k},\hat{\beta}_{2,\ell},\hat{\beta}_{2,k\ell})^{t}$,
respectively. Further, let
$\hat{\beta}_{2}^{\perp}=(\hat{\beta}_{2,\ell},\hat{\beta}_{2,k\ell})^{t}$ be
the estimated regression parameters of the covariates in model $M_{2}$, that
are additional to the covariates in model $M_{1}$. Finally, let
${\beta}_{1}^{0}=(\beta_{1,1}^{0},0)^{t}$, with
$\beta_{1,1}^{0}=\beta_{0,1}^{0}$ the parameter value in the true model
$M_{true}$.

######  Theorem 3.1 (Pairwise independence in the Cox PH model).

Let $(T_{1}=\tilde{T}_{1}\wedge
C_{1},\Delta_{1},{X}_{1}),\ldots,(T_{n}=\tilde{T}_{n}\wedge
C_{1},\Delta_{n},{X}_{n})$ be a random sample. Suppose that
$\tilde{T}_{i}|X_{i}$ is a random observation from the Cox PH model with
hazard function (17), then $\sqrt{n}(\hat{\beta}_{1}-{\beta}_{1}^{0})$ and
$\sqrt{n}\hat{\beta}_{2}^{\perp}$ are asymptotically independent.

The proof is given in Appendix D. Complete independence for the Cox PH model
can be shown as before, under the full null hypothesis and under the
assumptions of independent (and standardized) covariates. To control the FWER
in the two-stage testing procedures, the test-statistics in the two stages
need to be asymptotically completely independent, under the null hypothesis.
Like for the GLMs, this independence follows from the asymptotic (mixed)
pairwise independence of the maximum likelihood estimators of the regression
parameters and Slutky’s Lemma (see e.g., Van der Vaart (1998)).

##  4 Simulation studies and application

Suppose the data $(Y_{1},{X}_{1}),\ldots,(Y_{n},{X}_{n})$ of $n$ independent
individuals from a population of interest are available. For individual $i$
the variable $Y_{i}$ is the outcome of interest and
${X}_{i}=(x_{i,1},\ldots,x_{i,p})^{t}$ is a vector with a numerical
representation of the genotypes at $p$ biallelic markers. The genotypes can be
numerically represented in different ways depending on the genetic model that
is assumed (dominant, additive, recessive), but for simplicity we assume that
$x_{i,j}$ equals the number of minor alleles that individual $i$ has at marker
$j$; so taking the values 0, 1, 2. Our aim is to detect markers for which the
outcome of interest is associated with an interaction of markers.

In this section we present the results of simulation studies for the linear,
the Poisson and the Cox PH model. Type I error control of the whole two-stage
procedure is considered in Subsection 4.2, and in Subsection 4.3 power
comparisons for different settings (correlated and uncorrelated markers),
different first stage thresholds and for different values of the interaction
effects are made. In Subsection 4.4 we justify testing for marginal effects in
the first stage to detect interaction effects. We start with a description of
the simulation process in Subsection 4.1.

###  4.1 Simulation procedure

Data are simulated or obtained in two steps:

  1. 1.

obtaining marker data by simulation or using real life data,

  2. 2.

simulation of phenotype (outcome) data based on the marker data.

A single data set consists of marker and phenotype data of 2000 individuals
with each 3000 markers. Observations of different individuals are assumed to
be independent.

Marker data
The FWER and the power of the two stage testing procedure are studied in two
settings: the marker genotypes are either independent or correlated within
individuals.

Independent marker data. For the 2000 individuals with each 3000 markers the
number of minor alleles are simulated independently from a binomial
distribution with parameters $m=2$ (two alleles) and a probability that is
drawn from a uniform distribution on the interval $[0.1;0.5]$. For the power
simulations, two markers are considered to be causal. Their number of minor
alleles are drawn from the binomial distribution with parameters $m=2$ and
probability $0.2$. The minor allele frequencies of the causal markers are
fixed in order to diminish the variability in the data.

Correlated marker data. The correlated marker data are real life genotype data
to mimic a realistic correlation structure between the markers (Pecanka et al
(2019)). The data set contains the data of 2000 individuals with each 3000
markers. The linkage disequilibrium (LD) structure of the 3000 markers is
shown in Figure 1.

Figure 1: Heat map of the LD pattern for the 3000 markers based on data of
2000 individuals. LD is measured by the squared correlation coefficient.

A few quality checks have been performed to be sure that significant results
will not be the consequence of extremely low variability in the covariates or
the interaction term, or of collinearities in the second stage. In the data
sets with independent data hardly any marker failed to pass the quality check.

Phenotype data
Let $(x_{1},x_{2})$ be the pair of causal markers. For given values of main
and interaction effects $\beta_{1},\beta_{2}$ and $\beta_{3}$, the phenotype
(outcome) data are simulated from the models:

  * •

Linear model. The outcome value is simulated from the model
$Y=\beta_{1}x_{1}+\beta_{2}x_{2}+\beta_{3}x_{1}x_{2}+\epsilon$ by simulating
$\epsilon$ from the ${\mathcal{N}}(0,1)$ distribution.

  * •

Poisson model with offset. The offset $t$ is simulated from a uniform[1, 5]
distribution. Next, the outcome variable is simulated from a Poisson
distribution with mean
$\mu=t\;\exp(\beta_{1}x_{1}+\beta_{2}x_{2}+\beta_{3}x_{1}x_{2})$.

  * •

Cox PH model. The hazards function is assumed to be equal to
$\lambda(t)=\lambda_{0}(t)\exp(\beta_{1}x_{1}+\beta_{2}x_{2}+\beta_{3}x_{1}x_{2})$
with the baseline hazard $\lambda_{0}$ so that the event times follow a
Weibull distribution with shape and scale parameters equal to 2 and 1,
respectively, if $x_{1}=x_{2}=0$. The censoring times are simulated from a
uniform distribution at the interval $[z_{\\{0.70\\}};z_{\\{0.99\\}}]$ with
$z_{\\{x\\}}$ the $x$th quantile of the Weibull distribution with shape and
scale parameters equal to 2 and 1. These numbers were chosen so that
approximately 10% of the observations are censored under the full null
hypothesis that $\beta_{1}=\beta_{2}=\beta_{3}=0$.

Different values of the regression parameters $\beta_{1},\beta_{2}$ and
$\beta_{3}$ have been considered. In the simulation studies to study FWER
control, $\beta_{1}=\beta_{2}=\beta_{3}=0$.

###  4.2 Simulation results: FWER

For each model and multiple first stage thresholds, the FWER is studied for
independent and dependent marker genotypes. This is done in the following
steps: 1) if the markers are assumed to be independent, marker data are
simulated as described before (otherwise the correlated marker data are used).
2) phenotype (outcome) data are simulated for 2000 individuals under the full
null hypothesis of no main and no interaction effects
($\beta_{1}=\beta_{2}=\beta_{3}=0$). 3) the two-stage testing procedure is
applied. A type I error is made if at least one significant interaction is
found in the second stage. An interaction between two markers is significant
in the second stage if its p-value $<0.05/K_{1}$ with $K_{1}$ the number of
pairs in the second stage. 4) the steps 2 and 3 are repeated 5000 times. The
FWER is estimated by the fraction of data sets for which at least one
significant interaction was found in the second stage. This procedure is
performed for every model and several thresholds.

The results for the three models, for uncorrelated and correlated markers and
for the first stage thresholds equal to 0.05, 0.01 and 0.005 are presented in
Table 1. In case of uncorrelated markers, the FWERs are just below the
threshold of 0.05 for the three models and different first stage thresholds.
For the correlated markers, the FWERs are rather low; possibly the multiple
testing correction is too conservative. In Kooperberg et al (2008) simulation
studies had been performed for the logistic regression model. When studying
the type I error, they simulated the phenotype data from a model that includes
a main effect, but no interaction effect. The marker genotypes were simulated
from a first order Markov chain with a correlation between neighboring markers
equal to 0.7. They also found that the Bonferoni correction is conservative.
In a first order Markov chain, the correlation between the markers decreases
fast with their distance. That might be the reason that our results are more
conservative than their outcomes.

| uncorrelated markers | correlated markers
---|---|---
FST | Linear | Poisson | CoxPH | Linear | Poisson | CoxPH
0.05 | 0.0448 | 0.0416 | 0.0432 | 0.0076 | 0.0048 | 0.0064
0.01 | 0.0462 | 0.0492 | 0.0432 | 0.0114 | 0.0094 | 0.0120
0.005 | 0.0452 | 0.0364 | 0.0490 | 0.0078 | 0.0078 | 0.0084
Table 1: Type I error for different first stage thresholds (FST) for
uncorrelated and correlated markers.

###  4.3 Simulation results: Power

For this simulation study two causal markers are chosen. In the "uncorrelated
marker setting" the number of minor alleles at these causal markers are
sampled from a binomial distribution with parameters $m=2$ and probability
0.2. The number of minor alleles for the remaining markers are sampled as was
described before. For the correlated marker setting, two of the 3000 markers
are chosen such that they have a correlation of 0.3 and minor allele
frequencies equal to 0.48 and 0.25. These causal markers are used to simulate
the phenotype data from the three models of interest as described in
Subsection 4.1 with $x_{1}$ and $x_{2}$ the number of minor alleles at the two
causal markers. The interaction effect $\beta_{3}$ runs from 0 to 1 in small
steps of 0.05. Three settings for the main-effects are considered

  1. 1.

no main effects: $\beta_{1}=\beta_{2}=0$.

  2. 2.

two main effects: $\beta_{1}=\beta_{2}=0.5$.

  3. 3.

opposite main effects: $\beta_{1}=\beta_{2}=-0.5$.

The remaining simulation assumptions are as described in Subsections 4.1. For
each model, each value of the interaction effect ($\beta_{3}$) and for each
first stage threshold (FST $\in\\{1,0.1,0.05,0.01,0.005\\}$) new data are
simulated and the two stage testing procedure is performed. The power is
estimated as the fraction of times the causal pair of markers is significant
in the second stage after Bonferroni correction.

The results of the simulation studies for the uncorrelated markers in the
linear and the Cox model are presented in Figure 2. The results of the Poisson
regression model can be found in Appendix E. If there are no main effects
(only interaction effects, first column) or main effects in the same direction
as the interaction effect (second column) the two stage approach is clearly
more powerful than the approach in which all marker-marker combinations are
tested for interaction (i.e., FST=1). In general (except for one setting) the
power increases with a decreasing first stage threshold (FST). The lower the
threshold in the first stage, the less markers are passed on to the second
stage, and the less interaction terms are tested for association. This leads
to a less severe multiple testing correction. However, the lower the threshold
in the first stage, the lower the probability the causal markers are passed on
to the second stage. The power plots for opposite main effects show a drop in
power for large values of the interaction effect. This drop is not seen in the
classical approach (FST=1) and therefore, for large effects, the strategy of
testing all interaction pairs is more powerful than the two-stage testing
procedure. It seems that large interaction effects cancel out the marginal
association of the variables with the outcome, and prevent the markers to be
selected in the first stage.

Figure 2: Power as a function of the interaction effect for different first
stage threshold (FST), with independent markers. First row: linear regression
model; second row: Cox PH model. First column: no main effects; Second column:
both main effects equal 0.5; Third column: both main effects equal -0.5. (Note
the different scales on the x-axes.)

For the correlated marker data, the results for the linear regression model
and the Cox model are shown in Figure 3 and for the Poisson regression model
they can be found in Appendix E. The power gain of the two-stage testing
procedure is smaller in the correlated setting compared to the setting with
uncorrelated markers, but is still clearly visible, except in the setting in
which the marginal effects are opposite to the interaction effect. For first
stage thresholds smaller than 1, the power curves corresponding to the
different first stage threshold are almost equal.

A quality check has been performed in the second stage. Markers that are
strongly mutually correlated or with their interaction term ($r^{2}\geq 0.9$)
are not tested for an interaction effect and are left out from the analysis.
The check was performed to prevent unreliable results due to collinearity. The
two causal markers have a low correlation ($r^{2}=0.3$).

Figure 3: Power as a function of the interaction effect for different first
stage threshold (FST), with correlated markers. First row: linear regression
model; second row: Cox PH model. First column: no main effects; Second column:
both main effects equal 0.5; Third colum: both main effects equal -0.5. (Note
the different scales on the x-axes.)

###  4.4 Justification of using a marginal model for screening

In the first stage, markers are screened for interaction one by one by testing
for association between the marker and the outcome variable marginally. Even
if the marker is associated with the outcome via an interaction, and not
marginally, the test in this misspecified marginal model has power to detect
markers that are involved in an interaction. Mathematically this can be shown
in the linear model in which explicit expressions for the estimators exist,
but this is more complex to show in the logistic, Poisson and Cox PH models.
As an alternative we performed a simulation study. We simulated marker data
(the covariates $x_{1},x_{2}$) by sampling them independently from a binomial
distribution with parameters $m=2$ and probability 0.3 (to reflect the marker
data). Next we simulated outcome data for the different models that include
the interaction $x_{1}x_{2}$ only (no main effects). The simulation procedure
is as described before. Simulations are performed for interaction effects
$\beta$ running from -1 to 1 in small steps of 0.002. For every value of
$\beta$ marker and outcome data of 2000 individuals are sampled. For every
data set the Wald test-statistics is computed for testing for association
between $x_{1}$ and the outcome variable in a marginal regression model that
includes the variable $x_{1}$ only. In Figure 4 the value of the test-
statistic is plotted against the value of $\beta$ (the interaction-effect). In
all three plots it can be seen that the test-statistic is close to zero when
$\beta$ is close to zero and the value of the test-statistic deviates more
from zero (and the smaller the accompanying p-values) the more $\beta$ differs
from zero. This illustrates that the test in the first stage has power to
detect an interaction effect and justifies the choice of the test in the first
stage.

Figure 4: Scatter plot for the value of the Wald test-statistics of the main
effect in the marginal model against the value of the interaction effect in
the full model. Left: linear regression model, middle: Cox PH model, right:
Poisson regression model.

###  4.5 Data analysis

Non-muscle invasive bladder cancer (NMIBC) often recurs after initial
diagnosis. To investigate which markers are associated with the time from
initial diagnosis to recurrence, data from 1451 patients with NMIBC were
analyzed (Hof et al (2023)). For all patients, genotypes of 670.542 markers
are available for analysis. Of the 1451 patients, 591 had a recurrence during
their follow-up time. In this paper the aim is to identify marker-marker
interactions that are associated with the time to recurrence after initial
diagnosis with the proposed two-stage method. The significance threshold in
the first stage was chosen to be equal to 0.0001, later also 0.001 was
considered. When using 0.0001 none of the marker pairs was significant in the
second stage. For a first stage significance threshold of 0.001, the pair
rs4415228 with rs2567720 was the only combination that was significant after
Bonferroni correction, the corrected p-value was equal to 0.026. We have not
found any relevant biological meaning of these SNPs. As far as we know, these
SNPs have not been found in other bladder cancer studies.

##  5 Discussion

In the paper we described the methodology for a two stage testing procedure
with independent stages for detecting interactions in GLMs and the Cox
regression model. In the first stage all covariates are tested for association
with the outcome of interest in a model that includes the main effect of a
single covariate only. If the corresponding p-value is below a pre-specified
threshold, the covariate is passed on to the second stage. In the second stage
all pairs of covariates that survived the first stage are tested for
interaction.

The significance threshold in the first stage of the two-stage testing
procedure can be chosen freely. In the most extreme case the threshold is
taken equal to 1 for some covariates which implies that these covariates are
passed on to the second stage without being tested. In, for instance, in Dai
et al (2012) the researchers are interested in testing interaction between
genetic markers and environmental factors. They chose to screen the genetic
markers only in the first stage and to pass on the environmental factors to
the second stage directly; their threshold was taken equal to 1. For this
specific setting some of the proofs can be simplified.

In the simulation study we considered the setting with independent and with
correlated marker genotypes. For independent markers the type I errors were
close to the significance level of the test. In the correlated setting the
multiple testing correction was too conservative; the type I error was clearly
below the threshold. A less conservative multiple testing correction may lead
to more power. We leave this for a next project.

Although the two-stage testing strategy was proposed as a testing procedure
for detecting interactions, it may also detect main-effects via the
interaction term. By the choice of the test-statistics in the verification
stage, the testing procedure may not be very powerful for detecting main
effects. This is not seen as a weakness of the procedure as the focus is on
detecting interaction effects and main effects are by-catch.

In the simulation studies and the data analysis the calculations of the
p-values are based on the asymptotic normality of the test-statistic under the
null hypothesis. The normality approximations may not be very accurate for
extreme large and small (negative) values of the test-statistics. Since the
significance thresholds are usually very low, the low p-values in the second
stage can be recalculated with the Firth approximation to increase the
accuracy (see Firth (1993)). This will increases the calculation time, but
this does not need to be a computational burden as only a fraction of all
possible interactions are tested in the second stage. In the analysis of the
bladder cancer data, recalculation did not change the results.

As an application of the methodology, detection of interactions between
genetic markers or between genetic markers and environmental factors is given.
However, its applicability is much wider, it can be applied in any high
dimensional setting in which one is interested in detecting interaction
effects.

#### Acknowledgement

We would like to thank the Nijmegen Bladder Cancer Study project team of the
Department for IQ Health of the Radboudumc in Nijmegen for giving us access to
their database. This work made use of the Dutch national e-infrastructure with
the support of the SURF Cooperative using grant no. EINF-1968.

#### Software

An R package “twostageGWASsurvival” for doing the analysis have been developed
and is available at github: https://github.com/luc-van-
schijndel/twostageGWASsurvival. The package implements a two stage method for
the Cox model for searching for interaction effects within a high-dimensional
setting where the covariates outnumber the subjects.

## References

  * 1 Aalen O, Borgan O and Gjessing H.  Survival and Event History Analysis. Springer-Verlag, 2008.
  * 2 Andersen PK, Borgan O, Gill RD, Keiding N. Statistical Models Based on Counting Processes. Springer-Verlag, New York 1993.
  * 3 Dai JY, Kooperberg C, Leblanc M, Prentice RL  Two-stage testing procedures with independent filtering for genome-wide gene-environment interaction.  Biometrika 99(4): 929–944, 2012.
  * 4 Firth D,  Bias reduction of maximum likelihood estimates.  Biometrika 80, 27–38. 1993. doi: 10.1093/biomet/80.1.27
  * 5 Gauderman WJ, Thomas DC, Murray CE, Conti D, Li D, Lewinger J.  Efficient genome-wide association testing of gene-environment interaction in case-parent trios.  Am. J. Epidemiol. 172, 2010,
  * 6 Hof JP, Vermeulen SH, van der Heijden AG, Verhaegh GW, Dyrskjøt L, et al  An Association Study of Germline Variants in Bladder Cancer-Related Genes with the Prognosis of Non-Muscle Invasive Bladder Cancer.  Bladder Cancer, 9(1), 59-71, 2023
  * 7 Hsu L, Jiao S, Dai JY, Hutter C, Peters U, Kooperberg C.  Powerful Cocktail Methods for Detecting Genome-WideGene-Environment Interaction.  Genetic Epidemiology, 36: 183–194, 2012
  * 8 Johnsen PV, Riemer-Sørensen S, DeWan AT, Cahill ME, Langaas M.  A new method for exploring gene–gene and gene–environment interactions in GWAS with tree ensemble methods and SHAP values  BMC Bioinformatics, 22 (1), 2021.
  * 9 Kawaguchi ES, Li G, Lewinger JP, Gauderman WJ.  Two-step hypothesis testing to detect gene-environment interactions in a genome-wide scan with a survival endpoint.  Stat Med 41:1644–1657, 2022.
  * 10 Kawaguchi ES, Kim AE, Lewinger JP, Gauderman WJ.  Improved two-step testing of genome-wide gene–environment interactions  Genetic Epidemiology 47:152–166, 2023
  * 11 Kooperberg C, LeBlanc M  Increasing the power of identifying genexgene interactions in genome-wide association studies.  Genet Epidemiol, 32(3):255–263, 2008.
  * 12 Lewinger JP, Morrison JL, Thomas DC, Murcray CE, Conti DV, Li D, Gauderman WJ  Efficient two-step testing of gene-gene interactions in genomewide association studies.  Genet Epidemiol, 37(5):440–451, 2013.
  * 13 Murcray CE, Lewinger JP, Gauderman, WJ.  Gene–environment interaction in genome-wide associationstudies.  American Journal of Epidemiology, 169(2), 219–226. 2009. https://doi.org/10.1093/aje/kwn353.
  * 14 Pecanka J, Jonker MA, International Parkinson’s Disease Genomics Consortium (IPDGC), Bochdanovits Z, van der Vaart AW  A powerful and efficient two-stage method for detecting gene-to-gene interactions in GWAS.  Biostatistics 18(3):477–494, 2017
  * 15 Pecanka J, van der Vaart AW, Jonker MA.  Modeling association between multivariate correlated outcomes and high-dimensional sparse covariates: the adaptive SVS method.  Journal of Applied Statistics. 46 (5), 893-913, 2019.
  * 16 Wong K-C  Epistasis - Methods and Protocols. Chaper 6: Two-Stage Testing for Epistasis: Screening and Verification Pecanka J, Jonker MA  Methods in Molecular Biology 2212. Springer Protocols 2021\.
  * 17 Therneau TM, Grambsch PM.  Modeling survival data, extending the Cox model.  Springer-Verlag New York, Inc. 2000.
  * 18 van der Vaart AW.  Asymptotic Statistics.  Cambridge University Press, Cambridge, 1998.
  * 19 van der Vaart AW, and Wellner JA.  Weak convergence and Empirical processes. Springer-Verlag. New York. 1996.

Outline of the appendix:

  * -

Appendix A: Pairwise independence in GLMs

  * -

Appendix B: Complete independence in GLMs

  * -

Appendix C: FWER control

  * -

Appendix D: Pairwise independence in the Cox PH regression model

  * -

Appendix E: Power plots for the Poisson regression model

## Appendix A: Pairwise independence in GLMs

Theorem 2.1. Suppose that $Y_{1}|{{X}}_{1},\ldots,Y_{n}|{{X}}_{n}$ is a sample
from the generalized linear model with a canonical link function, $M_{true}$.
Let $\hat{{\theta}}_{1}$ and $\hat{{\theta}}_{2}$ denote the maximum
likelihood estimators of ${\theta}_{1}^{0}$ and ${{\theta}}_{2}^{0}$ in the
nested models $M_{1}$ and $M_{2}$, respectively. If model $M_{2}$ coincides
with the true model $M_{true}$ (i.e., $p_{2}=p_{0}$), then
$\sqrt{n}(\hat{\theta}_{1}-{\theta}_{1}^{0})$ and
$\sqrt{n}(\hat{\theta}_{2}^{\perp}-{\theta}_{2}^{0,\perp})$ are asymptotically
independent.

######  Proof .1.

The conditional density $f_{{\theta}_{0}}$ of $Y|{X}$ from the true model
$M_{true}$ is assumed to belong to a GLM with a canonical link function and,
therefore, can be written in the form:

$\displaystyle
f_{{\theta}_{0}}(Y|{X})=g(Y)\exp\big(V(Y)^{t}{Z}{\theta}_{0}-A({Z}{\theta}_{0})\big{missing}),$

where ${Z}$ is a function of the covariates ${X}$. The log-likelihood function
(divided by $n$) for the true model $M_{true}$ equals:

$\displaystyle\frac{1}{n}\ell({\theta}_{0})\;\propto\;\frac{1}{n}\sum_{i=1}^{n}\log\\{f_{{\theta}_{0}}(Y_{i}|{{X}}_{i})\big{\\}}\;=\;\frac{1}{n}\sum_{i=1}^{n}\big{\\{}V(Y_{i})^{t}{Z}_{i}{\theta}_{0}-A({Z}_{i}{\theta}_{0})\big{\\}}+\frac{1}{n}\sum_{i=1}^{n}\log\\{g(Y_{i})\big{\\}}.$

On the far right side, only the first term depends on the parameter of
interest ${\theta}_{0}$. Therefore, we can leave out the second term from the
calculations. Next, we consider the average (conditional) log-likelihood for
the submodel $M_{1}$:

$\displaystyle\frac{1}{n}\ell_{1}({\theta}_{1})\;\mathrel{\mathop{\mathchar
58\relax}}=\;\frac{1}{n}\ell(E_{1}{\theta}_{1})\;\propto\;\frac{1}{n}\sum_{i=1}^{n}\big{(}V(Y_{i})^{t}{Z}_{i}E_{1}{\theta}_{1}-A({Z}_{i}E_{1}{\theta}_{1})\big{)}.$
(18)

By the strong law of large numbers $\ell_{1}({\theta}_{1})/n$ converges
pointwise a.s. to its expectation:
$\mathbb{E}_{Y,{X}}\\{V(Y)^{t}{Z}E_{1}{\theta}_{1}-A({Z}E_{1}{\theta}_{1})\\}$.
Let ${\theta}_{1}^{0}\in\mathbb{R}^{q_{1}}$ be the maximizer of this limit. It
is known that the log-likelihood of an exponential family is concave,
therefore this maximizer is unique and a zero of its derivative with respect
to ${\theta}_{1}$. For $DA(.)$ defined as the gradient of $A(.)$ it holds
that:

$\displaystyle 0$
$\displaystyle=\;\mathbb{E}_{Y,{X}}\Big{\\{}V(Y)^{t}{Z}E_{1}-DA({Z}E_{1}{\theta}_{1}^{0}){Z}E_{1}\Big{\\}}$
$\displaystyle=\;\mathbb{E}_{{X}}\mathbb{E}_{Y|{X}}\Big{\\{}V(Y)^{t}{Z}E_{1}-DA({Z}E_{1}{\theta}_{1}^{0}){Z}E_{1}\Big{\\}}\;=\;\mathbb{E}_{{X}}\Big{\\{}DA({Z}{\theta}^{0}_{0}){Z}E_{1}-DA({Z}E_{1}{\theta}_{1}^{0}){Z}E_{1}\Big{\\}},$
(19)

with ${\theta}^{0}_{0}$ the true parameter value in model $M_{true}$. The last
equality follows from the exponential family identity:
$\mathbb{E}_{Y|{X}}(V(Y))=DA({Z}{\theta}^{0}_{0})^{t}$. Note that if
${\theta}^{0}_{0}=E_{1}{\theta}_{1}^{0}$, the right hand side is trivially
equal to zero.

By the argmax continuous mapping theorem and implications (see e.g., Van der
Vaart and Wellner (1996): Subsection 3.2) the maximum likelihood estimator for
$\theta_{1}$ in model $M_{1}$, the maximizer of the log-likelihood (18),
denoted as $\hat{\theta}_{1}$, converges in probability to ${\theta}_{1}^{0}$.
For any finite sample define $\hat{h}_{1}\mathrel{\mathop{\mathchar
58\relax}}=\hat{\theta}_{1}-{\theta}^{0}_{1}\in\mathbb{R}^{q_{1}}$ and
$\tilde{\ell}_{1}(\hat{h}_{1})\mathrel{\mathop{\mathchar
58\relax}}=\ell_{1}({\theta}_{1}^{0}+\hat{h}_{1})$. Then, by a Taylor
expansion of $A({Z}_{i}E_{1}({\theta}_{1}^{0}+\hat{h}_{1}))$ with respect to
$\hat{h}_{1}$ near zero:

$\displaystyle\frac{1}{n}\tilde{\ell}_{1}(\hat{h}_{1})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\Big{(}V(Y_{i})^{t}{Z}_{i}E_{1}({\theta}_{1}^{0}+\hat{h}_{1})-A({Z}_{i}E_{1}({\theta}_{1}^{0}+\hat{h}_{1}))\Big{)}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\Big{(}V(Y_{i})^{t}{Z}_{i}E_{1}\theta_{1}^{0}+V(Y_{i})^{t}{Z}_{i}E_{1}\hat{h}_{1}-A({Z}_{i}E_{1}\theta_{1}^{0})-DA({Z}_{i}E_{1}\theta_{1}^{0}){Z}_{i}E_{1}\hat{h}_{1}$
$\displaystyle~{}~{}~{}~{}\qquad-\frac{1}{2}\hat{h}_{1}^{t}E_{1}^{t}{Z}_{i}^{t}D^{2}A({Z}_{i}E_{1}\theta_{1}^{0}){Z}_{i}E_{1}\hat{h}_{1}\Big{)}+O_{p}\big{(}\|\hat{h}_{1}\|^{3}\big{)},$
(20)

where $DA$ and $D^{2}A$ are the first and second derivative of $A$ or the
gradient and the Hessian matrix.

The term $O_{p}(\|\hat{h}_{1}\|^{3})$ will converge to zero in probability.
The first and the third term within the summation in (20) do not depend on
$\hat{h}_{1}$ and their sum is written as $K_{n}$. By the law of large
numbers, the fifth term equals

$\displaystyle\frac{1}{2n}\sum_{i=1}^{n}\hat{h}_{1}^{t}E_{1}^{t}{Z}_{i}^{t}D^{2}A({Z}_{i}E_{1}\theta_{1}^{0}){Z}_{i}E_{1}\hat{h}_{1}=\tfrac{1}{2}\hat{h}_{1}^{t}E_{1}^{t}I_{1}E_{1}\hat{h}_{1}+o_{p}\big{(}\|\hat{h}_{1}\|^{2}\big{)},$

with $I_{1}\mathrel{\mathop{\mathchar
58\relax}}=\mathbb{E}_{X}{{Z}^{t}D^{2}A({Z}E_{1}{\theta}_{1}^{0}){Z}}\in\mathbb{R}^{q_{0}\times
q_{0}}$. The sum of the second and fourth term in (20) equals

$\frac{1}{n}\sum_{i=1}^{n}\Big{(}V(Y_{i})^{t}{Z}_{i}-DA({Z}_{i}E_{1}\theta_{1}^{0}){Z}_{i}\Big{)}E_{1}\hat{h}_{1}=\frac{1}{\sqrt{n}}{W}_{1}^{t}E_{1}\hat{h}_{1}+o_{p}\big{(}\|n^{-1/2}\hat{h}_{1}\|\big{)},$

by the central limit theorem, for ${W}_{1}\in\mathbb{R}^{q_{0}}$ normally
distributed.

The same derivation holds for submodel $M_{2}$, yielding the normally
distributed variable ${W}_{2}\in\mathbb{R}^{q_{0}}$. Jointly,

$\mathbb{R}^{1\times
2q_{0}}\ni{{W}_{1}\choose{W}_{2}}=\lim_{n\to\infty}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}{{Z}_{i}^{t}V(Y_{i})-{Z}_{i}^{t}DA({Z}_{i}E_{1}\theta_{1}^{0})^{t}\choose{Z}_{i}^{t}V(Y_{i})-{Z}_{i}^{t}DA(Z_{i}E_{2}\theta_{2}^{0})^{t}}$

has a normal distribution with mean zero by the equality in (19). Next we
compute its covariance-matrix. We start with the variance of $W_{1}$:

$\displaystyle\mathbb{E}_{Y,{X}}{{W}_{1}{W}_{1}^{t}}$
$\displaystyle=\mathbb{E}_{Y,{X}}\Big{\\{}\Big{(}{Z}^{t}V(Y)-{Z}^{t}DA(ZE_{1}\theta_{1}^{0})^{t}\Big{)}\Big{(}V(Y)^{t}{Z}-DA({Z}E_{1}\theta_{1}^{0})\Big{)}\Big{\\}}$
$\displaystyle=\mathbb{E}_{Y,{X}}\Big{\\{}\Big{(}\left({Z}^{t}V(Y)-{Z}^{t}DA({Z}\theta_{0}^{0})^{t}\right)-\left({Z}^{t}DA({Z}E_{1}\theta_{1}^{0})^{t}-{Z}^{t}DA(Z\theta_{0}^{0})^{t}\right)\Big{)}$
$\displaystyle\qquad\qquad\times\Big{(}{\left(V(Y)^{t}{Z}-DA({Z}\theta_{0}^{0}){Z}\right)-\left(DA({Z}E_{1}\theta_{1}^{0})Z-DA({Z}\theta_{0}^{0})Z\right)}\Big{)}\Big{\\}}$
$\displaystyle=\mathbb{E}_{Y,{X}}\Big{\\{}\Big{(}{Z}^{t}V(Y)-{Z}^{t}DA({Z}\theta_{0}^{0})^{t}\Big{)}\Big{(}V(Y)^{t}{Z}-DA({Z}\theta_{0}^{0})Z\Big{)}\Big{\\}}$
$\displaystyle\qquad+\mathbb{E}_{Y.{X}}\Big{\\{}\Big{(}{Z}^{t}DA({Z}E_{1}\theta_{1}^{0})^{t}-Z^{t}DA({Z}\theta_{0}^{0})^{t}\Big{)}\Big{(}DA(ZE_{1}\theta_{1}^{0}){Z}-DA({Z}\theta_{0}^{0})Z\Big{)}\Big{\\}}$
$\displaystyle\qquad-\mathbb{E}_{Y,{X}}\Big{\\{}\Big{(}{Z}^{t}V(Y)-{Z}^{t}DA({Z}\theta_{0}^{0})^{t}\Big{)}\Big{(}DA({Z}E_{1}\theta_{1}^{0}){Z}-DA({Z}\theta_{0}^{0})\Big{)}\Big{\\}}$
$\displaystyle\qquad-\mathbb{E}_{Y,{X}}\Big{\\{}\Big{(}{Z}^{t}DA({Z}E_{1}\theta_{1}^{0})^{t}-{Z}^{t}DA({Z}\theta_{0}^{0})^{t}\Big{)}\Big{(}V(Y)^{t}{Z}-DA({Z}\theta_{0}^{0})\Big{)}\Big{\\}}.$

The first of the four terms can be further simplified, by applying the
exponential family-identity,
$\mbox{Cov}_{Y|{X}}(V(Y))=D^{2}A(Z\theta_{0}^{0})$:

$\displaystyle\
\mathbb{E}_{Y,{X}}\Big{\\{}\Big{(}Z^{t}V(Y)-{Z}^{t}DA(Z\theta_{0}^{0})^{t}\Big{)}\Big{(}V(Y)^{t}Z-DA(Z\theta_{0}^{0}){Z}\Big{)}\Big{\\}}$
$\displaystyle\qquad=\
\mathbb{E}_{X}{Z^{t}\mathbb{E}_{Y|{X}}\Big{\\{}\Big{(}V(Y)-DA({Z}\theta_{0}^{0})^{t}\Big{)}\Big{(}V(Y)-DA({Z}\theta_{0}^{0})^{t}}\Big{)}^{t}\Big{\\}}$
$\displaystyle\qquad=\
\mathbb{E}_{X}\Big{\\{}{Z^{t}\mbox{Cov}_{Y|X}\big{(}V(Y)\big{)}}\Big{\\}}$
$\displaystyle\qquad=\
\mathbb{E}_{X}\Big{\\{}Z^{t}D^{2}A(Z\theta_{0}^{0}){Z}\Big{\\}}\;=\mathrel{\mathop{\mathchar
58\relax}}\ I\in\mathbb{R}^{q_{0}\times q_{0}}.$ (21)

The second term can not be simplified, we define it as:

$\displaystyle C_{11}$ $\displaystyle\mathrel{\mathop{\mathchar
58\relax}}=\mathbb{E}_{X}\Big{\\{}\Big{(}Z^{t}DA(ZE_{1}\theta_{1}^{0})^{t}-Z^{t}DA(Z\theta_{0}^{0})^{t}\Big{)}\Big{(}DA(ZE_{1}\theta_{1}^{0})Z-DA(Z\theta_{0}^{0})Z\Big{)}\Big{\\}}\in\mathbb{R}^{q_{0}\times
q_{0}}.$

The last two terms equal zero. This can be seen as follows:

$\displaystyle\
\mathbb{E}_{Y,X}\Big{\\{}\Big{(}Z^{t}V(Y)-Z^{t}DA(Z\theta_{0}^{0})^{t}\Big{)}\Big{(}DA(ZE_{1}\theta_{1}^{0})Z-DA(Z\theta_{0}^{0})Z\Big{)}\Big{\\}}$
$\displaystyle\qquad=\
\mathbb{E}_{X}\mathbb{E}_{Y|X}\Big{\\{}\Big{(}Z^{t}V(Y)-Z^{t}DA(Z\theta_{0}^{0})^{t}\Big{)}\Big{(}DA(ZE_{1}\theta_{1}^{0})Z-DA(Z\theta_{0}^{0})Z\Big{)}\Big{\\}}$
$\displaystyle\qquad=\
\mathbb{E}_{X}\Big{\\{}\Big{(}Z^{t}\mathbb{E}_{Y|X}\big{(}V(Y)\big{)}-Z^{t}DA(Z\theta_{0}^{0})^{t}\Big{)}\Big{(}DA(ZE_{1}\theta_{1}^{0})Z-DA(Z\theta_{0}^{0})Z\Big{)}\Big{\\}}\;=\;0,$

since $\mathbb{E}_{Y|X}\\{V(Y)\\}=DA(Z\theta_{0}^{0})^{t}$. For the fourth
term similar reasoning holds. Summarized

$\displaystyle\mathbb{E}_{Y,X}\big{\\{}W_{1}W_{1}^{t}\big{\\}}=I+C_{11}.$

Similarly, we find $\mathbb{E}_{X,Y}\\{W_{2}W_{2}^{t}\\}=I+C_{22}$, with
$C_{22}$ defined as

$\displaystyle C_{22}$ $\displaystyle\mathrel{\mathop{\mathchar
58\relax}}=\mathbb{E}_{X}\Big{\\{}\Big{(}Z^{t}DA(ZE_{2}\theta_{2}^{0})^{t}-Z^{t}DA(Z\theta_{0}^{0})^{t}\Big{)}\Big{(}DA(ZE_{2}\theta_{2}^{0})Z-DA(Z\theta_{0}^{0})Z\Big{)}\Big{\\}}\in\mathbb{R}^{q_{0}\times
q_{0}}.$

The covariance of $W_{1}$ and $W_{2}$:

$\displaystyle\mathbb{E}_{Y,X}{W_{1}W_{2}^{t}}$
$\displaystyle=\mathbb{E}_{Y,X}\Big{\\{}\Big{(}Z^{t}V(Y)-Z^{t}DA(ZE_{1}\theta_{1}^{0})^{t}\Big{)}\Big{(}V(Y)^{t}Z-DA(ZE_{2}\theta_{2}^{0})Z\Big{)}\Big{\\}}$
$\displaystyle=\mathbb{E}_{Y,X}\Big{\\{}\Big{(}\big{(}Z^{t}V(Y)-Z^{t}DA(Z\theta_{0}^{0})^{t}\big{)}-\big{(}Z^{t}DA(ZE_{1}\theta_{1}^{0})^{t}-Z^{t}DA(Z\theta_{0}^{0})^{t}\big{)}\Big{)}$
$\displaystyle\qquad\times\Big{(}\big{(}V(Y)^{t}Z-DA(Z\theta_{0}^{0})Z\big{)}-\big{(}DA(ZE_{2}\theta_{2}^{0})Z-DA(Z\theta_{0}^{0})Z\big{)}\Big{)}\Big{\\}}$
$\displaystyle=\mathbb{E}_{X}\Big{\\{}Z^{t}D^{2}A(Z\theta_{0}^{0})Z\Big{\\}}\;+\;\mathbb{E}_{X}\Big{\\{}\Big{(}Z^{t}DA(ZE_{1}\theta_{1}^{0})^{t}-Z^{t}DA(Z\theta_{0}^{0})^{t}\Big{)}\Big{(}DA(ZE_{2}\theta_{2}^{0})Z-DA(Z\theta_{0}^{0})Z\Big{)}\Big{\\}}$
$\displaystyle=I+C_{12}$

with (by similar reasoning)

$\displaystyle C_{12}\mathrel{\mathop{\mathchar
58\relax}}=\mathbb{E}_{X}\Big{\\{}\Big{(}Z^{t}DA(ZE_{1}\theta_{1}^{0})^{t}-Z^{t}DA(Z\theta_{0}^{0})^{t}\Big{)}\Big{(}DA(ZE_{2}\theta_{2}^{0})Z-DA(Z\theta_{0}^{0})Z\Big{)}\Big{\\}}\in\mathbb{R}^{q_{0}\times
q_{0}},$

and $C_{12}={C_{12}}^{t}$. Now, due to the combined vector
$(W_{1}^{t},W_{2}^{t})^{t}$ having expectation zero, we find its covariance
matrix:

$\mbox{Cov}_{Y,X}{W_{1}\choose
W_{2}}=\begin{pmatrix}I+C_{11}&~{}~{}~{}~{}I+C_{21}\\\
I+C_{12}&~{}~{}~{}~{}I+C_{22}\end{pmatrix}\in\mathbb{R}^{2q_{0}\times
2q_{0}}.$

Note that if the submodel $M_{2}$ coincides with model $M_{true}$, then
$E_{2}\theta_{2}^{0}=\theta_{0}^{0}$ and
$C_{21}=C_{12}=C_{22}=\mathbf{0}\in\mathbb{R}^{q_{0}\times q_{0}}$.

Summarized

$\displaystyle\frac{1}{n}\tilde{\ell}(\hat{h}_{1})=K_{n}+\frac{1}{\sqrt{n}}W_{1}^{t}E_{1}\hat{h}_{1}-\frac{1}{2}\hat{h}_{1}^{t}E_{1}^{t}I_{1}E_{1}\hat{h}_{1}+o_{p}\big{(}\|n^{-1/2}\hat{h}_{1}\|\big{)}+o_{p}\big{(}\|\hat{h}_{1}\|^{2}\big{)},$
(22)

In a similar way $n^{-1}\tilde{\ell}(\hat{h}_{2})$ can be written as a sum of
asymptotic terms, with $I_{2}\mathrel{\mathop{\mathchar
58\relax}}=\mathbb{E}_{X}{Z^{t}D^{2}A(ZE_{2}\theta_{2}^{0})Z}\in\mathbb{R}^{q_{0}\times
q_{0}}$. Maximizing $n^{-1}\tilde{\ell}(h_{1})$ with respect to $h_{1}$ to
obtain an expression for $\hat{h}_{1}$ gives for $\hat{h}_{1}$, and in a
similar way for $\hat{h}_{2}$:

$\displaystyle\hat{h}_{1}$
$\displaystyle=\frac{1}{\sqrt{n}}\left(E_{1}^{t}I_{1}E_{1}\right)^{-1}E_{1}^{t}W_{1}\in\mathbb{R}^{q_{1}}$
$\displaystyle\hat{h}_{2}$
$\displaystyle=\frac{1}{\sqrt{n}}\left(E_{2}^{t}I_{2}E_{2}\right)^{-1}E_{2}^{t}W_{2}\in\mathbb{R}^{q_{2}}.$

In order to analyze the asymptotic covariance of $\hat{h}_{1}$ and
$\hat{h}_{2}$ orthogonal to the submodel $M_{1}$, this is found by
$(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2})\hat{h}_{2}\in\mathbb{R}^{q_{2}}$,
where $Id_{q_{2}}$ denotes the $q_{2}$-dimensional identity matrix.
Ultimately, the $q_{1}\times q_{2}$-covariance becomes:

$\displaystyle\mbox{Cov}_{Y,X}\Big{(}$
$\displaystyle\hat{h}_{1},\left(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2}\right)\hat{h}_{2}\Big{)}$
$\displaystyle=\frac{1}{n}\left(E_{1}^{t}I_{1}E_{1}\right)^{-1}E_{1}^{t}\;\mbox{Cov}_{Y,X}(W_{1},W_{2})E_{2}\left(E_{2}^{t}I_{2}E_{2}\right)^{-1}\left(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2}\right)$
$\displaystyle=\frac{1}{n}\left(E_{1}^{t}I_{1}E_{1}\right)^{-1}E_{1}^{t}\left(I+C_{12}\right)E_{2}\left(E_{2}^{t}I_{2}E_{2}\right)^{-1}\left(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2}\right).$
(23)

To complete the proof, we set $M_{2}=M_{true}$, so
$\theta_{0}^{0}=E_{2}\theta_{2}^{0}$ which in turn implies $I_{2}=I$ and
$C_{21}=\mathbf{0}$. By the definitions of $E_{1}$ and $E_{2}$ as projection
matrices: Furthermore:

$\displaystyle E_{1}^{t}$ $\displaystyle=E_{1}^{t}E_{2}E_{2}^{t}$ (24)
$\displaystyle E_{1}^{t}$ $\displaystyle=E_{1}^{t}E_{1}E_{1}^{t},$ (25)

where the equality in (24) follows from the fact that $M_{1}$ is nested within
$M_{2}$. Apply these to find:

$\displaystyle\mbox{Cov}_{Y,X}\big{(}\sqrt{n}(\hat{\theta}_{1}-{\theta}_{1}^{0}),\sqrt{n}(\hat{\theta}_{2}-{\theta}_{2}^{0})^{\perp}\big{)}$
$\displaystyle=n\;\mbox{Cov}_{Y,X}\big{(}\hat{h}_{1},\left(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2}\right)\hat{h}_{2}\big{)}$
$\displaystyle=\left(E_{1}^{t}I_{1}E_{1}\right)^{-1}E_{1}^{t}E_{2}\left(E_{2}^{t}\left(I_{2}+\mathbf{0}\right)E_{2}\right)\left(E_{2}^{t}I_{2}E_{2}\right)^{-1}\left(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2}\right)$
$\displaystyle=\left(E_{1}^{t}I_{1}E_{1}\right)^{-1}E_{1}^{t}E_{2}\left(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2}\right)$
$\displaystyle=\left(E_{1}^{t}I_{1}E_{1}\right)^{-1}\left(E_{1}^{t}E_{2}-E_{1}^{t}E_{2}E_{2}^{t}E_{1}E_{1}^{t}E_{2}\right)$
$\displaystyle=\left(E_{1}^{t}I_{1}E_{1}\right)^{-1}\left(E_{1}^{t}E_{2}-E_{1}^{t}E_{1}E_{1}^{t}E_{2}\right)$
$\displaystyle=\left(E_{1}^{t}I_{1}E_{1}\right)^{-1}\left(E_{1}^{t}E_{2}-E_{1}^{t}E_{2}\right)$
$\displaystyle=\mathbf{0},$ (26)

where the second and fifth equalities follow from (24), and the sixth equality
follows from (25). By the joint normality of $\sqrt{n}\;\hat{h}_{1}$ and
$\sqrt{n}\;\hat{h}_{2}$, a zero correlation implies the required independence
as stated in the theorem.

Theorem 2.2 (Pairwise independence in the linear regression model). Suppose
that $Y_{1}|{{X}}_{1},\ldots,Y_{n}|{{X}}_{n}$ is a sample from the linear
regression model $M_{2}$ in (9) with the identical link function and unknown
variance $\sigma^{2}$. Then $\hat{\beta}_{1}$ and $\hat{\beta}_{2}^{\perp}$
are independent.

######  Proof .2.

The pair $(\hat{\beta}_{1},\hat{\beta}_{2})$ follows a multivariate normal
distribution as $\hat{\beta}_{1}$ and $\hat{\beta}_{2}$ are linear
combinations of $Y_{1},\ldots,Y_{n}$. To proof independence between
$\hat{\beta}_{1}$ and $\hat{\beta}_{2}^{\perp}$, we have to show that their
covariance equals zero:
$\mbox{Cov}(\hat{\beta}_{1},\hat{\beta}_{2}^{\perp})=0$. To prove this, we
compute the covariance $\mbox{Cov}(\hat{\beta}_{1},\hat{\beta}_{2})$ and show
that the elements of the covariance matrix that represent the covariance
between $\hat{\beta}_{1}$ and $\hat{\beta}_{2}^{\perp}$ equal zero.

Define ${Y}=(Y_{1},\ldots,Y_{n})^{t}$ and ${Z}_{1}$ and ${Z}_{2}$ as the
$p_{1}\times n$ and $p_{2}\times n$ matrices with the covariates of the $n$
observations in the models $M_{1}$ and $M_{2}$ in (4) and (5). Then,
$\hat{\beta}_{1}=({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Y}$ and
$\hat{\beta}_{2}=({Z}_{2}{Z}_{2}^{t})^{-1}{Z}_{2}{Y}$

$\displaystyle\mbox{Cov}\big{(}\hat{\beta}_{1},\hat{\beta}_{2}\big{)}$
$\displaystyle=\mbox{Cov}\big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Y},({Z}_{2}{Z}_{2}^{t})^{-1}{Z}_{2}{Y}\big{)}$
$\displaystyle=\mathbb{E}_{{Z},{Y}}\Big{(}({Z}_{1}{Z}_{1}^{T})^{-1}{Z}_{1}{Y}{Y}^{t}{Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1}\Big{)}-\mathbb{E}_{{Z},{Y}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Y}\Big{)}\mathbb{E}_{{Z},{Y}}\Big{(}{Y}^{t}{Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1}\Big{)}$

with

$\displaystyle\mathbb{E}_{{Z},{Y}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Y}\Big{)}\;$
$\displaystyle=\mathbb{E}_{{Z}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}\Big{)}\mathbb{E}_{{Y}|{Z}}{Y}\;=\;\mathbb{E}_{{Z}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Z}_{2}^{t}{\beta}_{2}^{0}\Big{)}$
$\displaystyle\mathbb{E}_{{Z},{Y}}\Big{(}{Y}^{t}{Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1}\Big{)}\;$
$\displaystyle=\mathbb{E}_{{Z}}\mathbb{E}_{{Y}|{Z}}\Big{(}{Y}^{t}{Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1}\Big{)}\;=\;\mathbb{E}_{{Z}}\Big{(}({\beta}_{2}^{0})^{t}{Z}_{2}{Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1}\Big{)}\;=\;({\beta}_{2}^{0})^{t}$
$\displaystyle\mathbb{E}_{{Z},{Y}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Y}\Big{)}\mathbb{E}_{{Z},{Y}}\Big{(}{Y}^{t}{Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1}\Big{)}\;$
$\displaystyle=\;\mathbb{E}_{{Z}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Z}_{2}^{t}{\beta}_{2}^{0}({\beta}_{2}^{0})^{t}\Big{)}$

and

$\displaystyle\mathbb{E}_{{Z},{Y}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Y}{Y}^{t}{Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1})\Big{)}$
$\displaystyle=\;\mathbb{E}_{{Z}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}\Big{)}\mathbb{E}_{{Y}|{Z}}\Big{(}({Y}{Y}^{t}){Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1}\Big{)}$
$\displaystyle=\;\mathbb{E}_{{Z}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}\big{\\{}\sigma^{2}I_{n}+({Z}_{2}^{t}{\beta}_{2}^{0})({Z}_{2}^{t}{\beta}_{2}^{0})^{t}\big{\\}}\;{Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1}\Big{)}$
$\displaystyle=\;\sigma^{2}\mathbb{E}_{{Z}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1}\Big{)}\;+\;\mathbb{E}_{{Z}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Z}_{2}^{t}{\beta}_{2}^{0}({Z}_{2}^{t}{\beta}_{2}^{0})^{t}{Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1}\Big{)}$
$\displaystyle=\;\sigma^{2}\mathbb{E}_{{Z}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1}\Big{)}\;+\;\mathbb{E}_{{Z}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Z}_{2}^{t}{\beta}_{2}^{0}({\beta}_{2}^{0})^{t}\Big{)}$

That gives

$\displaystyle\mbox{Cov}\big{(}\hat{\beta}_{1},\hat{\beta}_{2}\big{)}$
$\displaystyle=\sigma^{2}\mathbb{E}_{{Z}}\Big{(}({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1}\Big{)}$

Since model $M_{1}$ is nested in model $M_{2}$, we can write
${Z}_{2}^{t}=({Z}_{1}^{t}|{Z}_{3}^{t})$ with ${Z}_{3}$ the covariate rows
which are in ${Z}_{2}$ but not in ${Z}_{1}$. This gives
${Z}_{1}{Z}_{2}^{t}={Z}_{1}({Z}_{1}^{t}|{Z}_{3}^{t})=({Z}_{1}{Z}_{1}^{t}|{Z}_{1}{Z}_{3}^{t})$.
And

$\displaystyle({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Z}_{2}^{t}({Z}_{2}{Z}_{2}^{t})^{-1}\;=\;({Z}_{1}{Z}_{1}^{t})^{-1}({Z}_{1}{Z}_{1}^{t}|{Z}_{1}{Z}_{3}^{t})({Z}_{2}{Z}_{2}^{t})^{-1}\;=\;(I_{p_{1}}\;|\;({Z}_{1}{Z}_{1}^{t})^{-1}{Z}_{1}{Z}_{3}^{t})({Z}_{2}{Z}_{2}^{t})^{-1},$

with $I_{p_{1}}$ the $p_{1}\times p_{1}$ identity matrix. Write

$\displaystyle{Z}_{2}{Z}_{2}^{t}\;=\;\begin{pmatrix}{Z}_{1}{Z}_{1}^{t}&~{}~{}{Z}_{1}{Z}_{3}^{t}\\\
{Z}_{3}{Z}_{1}^{t}&~{}~{}{Z}_{3}{Z}_{3}^{t}\end{pmatrix}$

Then, to shorten notation in the next display, let
${W}_{11}\mathrel{\mathop{\mathchar 58\relax}}={Z}_{1}{Z}_{1}^{t}$,
${W}_{33}\mathrel{\mathop{\mathchar 58\relax}}={Z}_{3}{Z}_{3}^{t}$ and
${W}_{13}\mathrel{\mathop{\mathchar 58\relax}}={Z}_{1}{Z}_{3}^{t}$

$\displaystyle\begin{pmatrix}{Z}_{1}{Z}_{1}^{t}&~{}~{}{Z}_{1}{Z}_{3}^{t}\\\
{Z}_{3}{Z}_{1}^{t}&~{}~{}{Z}_{3}{Z}_{3}^{t}\end{pmatrix}^{-1}$
$\displaystyle=\;\begin{pmatrix}{W}_{11}&~{}~{}{W}_{13}\\\
{W}_{13}^{t}&~{}~{}{W}_{33}\end{pmatrix}^{-1}$
$\displaystyle=\begin{pmatrix}{W}_{11}^{-1}+{W}_{11}^{-1}{W}_{13}({W}_{33}-{W}_{13}^{t}{W}_{11}^{-1}{W}_{13})^{-1}{W}_{13}^{t}{W}_{11}^{-1}&~{}~{}~{}-{W}_{11}^{-1}{W}_{13}({W}_{33}-{W}_{13}^{t}{W}_{11}^{-1}{W}_{13})^{-1}\\\
-({W}_{33}-{W}_{13}^{t}{W}_{11}^{-1}{W}_{13})^{-1}{W}_{13}^{t}{W}_{11}^{-1}&({W}_{33}-{W}_{13}^{t}{W}_{11}^{-1}{W}_{13})^{-1}\end{pmatrix}$

We need to show that the last block column of the matrix

$\displaystyle(I_{p_{1}}\;|\;{W}_{11}^{-1}{W}_{13})$
$\displaystyle\begin{pmatrix}{W}_{11}&~{}~{}{W}_{13}\\\
{W}_{13}^{t}&~{}~{}{W}_{33}\end{pmatrix}^{-1}$

equals zero. That block equals

$\displaystyle-{W}_{11}^{-1}{W}_{13}({W}_{33}-{W}_{13}^{t}{W}_{11}^{-1}{W}_{13})^{-1}+{W}_{11}^{-1}{W}_{13}({W}_{33}-{W}_{13}^{t}{W}_{11}^{-1}{W}_{13})^{-1}={\bf
0}$

## Appendix B: Mixed pairwise independence in GLMs

In this Appendix we will prove Theorem 2.3. If in Theorem 2.3, $j=k$ or
$j=\ell$ the stated independence is pairwise independence and follows from
Result 2.1. Here we assume that $j\neq k$ and $j\neq\ell$. First, we give the
proof in a more general setting. After the proof, the theorem is translated to
the setting we are interested in. Define the models

$\displaystyle M_{1}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{X},{\beta}_{1},{\nu}_{1})\big{)}$
$\displaystyle=\beta_{1,1}x_{1}+\beta_{1,2}x_{2}+\ldots+\beta_{1,p_{1}}x_{p_{1}}$
(27) $\displaystyle M_{2}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{X},{\beta}_{2},{\nu}_{2})\big{)}$
$\displaystyle=\beta_{2,1}x_{1}+\beta_{2,2}x_{2}+\hskip
56.9055pt\beta_{2,p_{1}+1}x_{p_{1}+1}+\ldots+\beta_{2,p_{1}+p_{2}}x_{p_{1}+p_{2}}$
(28) $\displaystyle M_{3}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{X},{\beta}_{2},{\nu}_{3})\big{)}$
$\displaystyle=\beta_{3,1}x_{1}+\beta_{3,2}x_{2}+\ldots+\beta_{3,p_{1}}x_{p_{1}}+\beta_{3,p_{1}+1}x_{p_{1}+1}+\ldots+\beta_{3,p_{1}+p_{2}}x_{p_{1}+p_{2}}$
(29) $\displaystyle M_{true}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{X},{\beta}_{0},{\nu}_{0})\big{)}$
$\displaystyle=\beta_{0,1}x_{1}+\beta_{0,2}x_{2}$ (30)

where the models $M_{1},M_{2}$ and $M_{true}$ are defined as before. The
models $M_{1},M_{2}$ and $M_{true}$ are nested in model $M_{3}$. Model
$M_{true}$ is the true model with only the intercept and the “fixed”
covariates. Define the vectors $\hat{\theta}_{i}$ as the maximum likelihood
estimator of the full parameter vector in model $M_{i}$, and ${\theta}_{0}$ as
the true parameter vector (so the regression parameters and the nuisance
parameters). The dimension of the parameter spaces are denoted as
$q_{1},q_{2},q_{3}$ and $q_{0}$.

The parameter spaces for the models $M_{1}$ and $M_{2}$ are nested in model
$M_{3}$. To project a parameter from the parameter space for $M_{1}$ or
$M_{2}$ to the parameter space for model $M_{3}$ we define the projection
matrices (as before)

$\displaystyle E_{1}$
$\displaystyle=\left(e_{1}|\ldots|e_{q_{1}}\right)\in\mathbb{R}^{q_{3}\times
q_{1}}$ $\displaystyle E_{2}$
$\displaystyle=\left(e_{1}|\ldots|e_{q_{2}}\right)\in\mathbb{R}^{q_{3}\times
q_{2}},$

in such a way that $\\{e_{1},\ldots,e_{q_{1}}\\}$ forms an orthonormal basis
for the parameter space for $M_{1}$ and similarly for $E_{2}$. The matrices
represent linear functions from the parameter spaces for $M_{1}$ and $M_{2}$
to $M_{3}$.

######  Lemma .3 (Mixed pairwise independence in GLMs).

Suppose that $Y_{1}|{{X}}_{1},\ldots,Y_{n}|{{X}}_{n}$ is a sample from the
generalized linear model $M_{true}$ in (30) with a canonical link function.
Let $\hat{{\theta}}_{1}$ and $\hat{{\theta}}_{2}$ be the maximum likelihood
estimators of ${\theta}_{1}$ and ${{\theta}}_{2}$ in the models $M_{1}$ and
$M_{2}$ in (27) and (28), respectively. If
$E_{2}E_{2}^{t}\mathbb{E}_{X}{Z^{t}\mbox{Cov}_{Y|X}(V(Y))Z}E_{2}=\mathbb{E}_{X}{Z^{t}\mbox{Cov}_{Y|X}(V(Y))Z}E_{2}$,
where $Z$ and $V(Y)$ are defined as in the expression of the exponential
family (7), then $\sqrt{n}(\hat{\theta}_{1}-{\theta}_{1}^{0})$ and
$\sqrt{n}(\hat{\theta}_{2}^{0}-{\theta}_{2}^{0})^{\perp}$ are asymptotically
independent.

######  Proof .4.

The statement in this lemma can be proven along the same lines as the proof of
Theorem 2.1 up to the calculation of the covariance in (.1) at the end of the
proof. In that calculation the equality $E_{1}^{t}=E_{1}^{t}E_{2}E_{2}^{t}$ in
(24) is applied, but does not hold in the current setting as the model $M_{1}$
is not nested in $M_{2}$. However, if
$E_{2}E_{2}^{t}\mathbb{E}_{X}{Z^{t}\mbox{Cov}_{Y|X}(V(Y))Z}E_{2}=\mathbb{E}_{X}{Z^{t}\mbox{Cov}_{Y|X}(V(Y))Z}E_{2}$,
the covariance can be proven to be equal to zero in the same way. This is
explained below. Remind that $I=\mathbb{E}_{X}{Z^{t}\mbox{Cov}_{Y|X}(V(Y))Z}$.
Under the null hypothesis that there are no main and no interaction effects,
the true model coincides with the models $M_{1}$ and $M_{2}$. Then, the
covariance of interest is given by

$\displaystyle\mbox{Cov}_{Y,X}\big{(}\sqrt{n}(\hat{\theta}_{1}-{\theta}_{1}^{0}),\sqrt{n}(\hat{\theta}_{2}-{\theta}_{2}^{0})^{\perp}\big{)}$
$\displaystyle=n\;\mbox{Cov}_{Y,X}\big{(}\hat{h}_{1},\left(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2}\right)\hat{h}_{2}\big{)}$
$\displaystyle=\left(E_{1}^{t}IE_{1}\right)^{-1}E_{1}^{t}IE_{2}\left(E_{2}^{t}IE_{2}\right)^{-1}\left(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2}\right)$

Using the assumption that $IE_{2}=E_{2}E_{2}^{t}IE_{2}$ yields

$\displaystyle\left(E_{1}^{t}IE_{1}\right)^{-1}E_{1}^{t}IE_{2}\left(E_{2}^{t}IE_{2}\right)^{-1}\left(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2}\right)$
$\displaystyle=\left(E_{1}^{t}IE_{1}\right)^{-1}E_{1}^{t}E_{2}E_{2}^{t}IE_{2}\left(E_{2}^{t}IE_{2}\right)^{-1}\left(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2}\right)$
$\displaystyle=\left(E_{1}^{t}IE_{1}\right)^{-1}E_{1}^{t}E_{2}\left(Id_{q_{2}}-E_{2}^{t}E_{1}E_{1}^{t}E_{2}\right)$
$\displaystyle=\left(E_{1}^{t}IE_{1}\right)^{-1}E_{1}^{t}E_{2}-\left(E_{1}^{t}IE_{1}\right)^{-1}E_{1}^{t}E_{2}E_{2}^{t}E_{1}E_{1}^{t}E_{2}$
$\displaystyle=\left(E_{1}^{t}IE_{1}\right)^{-1}E_{1}^{t}E_{2}-\left(E_{1}^{t}IE_{1}\right)^{-1}E_{1}^{t}E_{1}E_{1}^{t}E_{2}$
$\displaystyle=\left(E_{1}^{t}IE_{1}\right)^{-1}E_{1}^{t}E_{2}-\left(E_{1}^{t}IE_{1}\right)^{-1}E_{1}^{t}E_{2}$
$\displaystyle=\mathbf{0}$

The lemma above is proven under the assumptions that
$E_{2}E_{2}^{t}\mathbb{E}_{X}{Z^{t}\mbox{Cov}_{Y|X}(V(Y))Z}E_{2}=\mathbb{E}_{X}{Z^{t}\mbox{Cov}_{Y|X}(V(Y))Z}E_{2}$.
By straightforward calculations, it can be seen that for the linear, the
Poisson, and the logistic regression model this assumption is satisfied if the
covariates that are included in the two stage testing procedure (so not
necessarily the fixed covariates) are mutually independent and independent of
the fixed covariates, and the covariates are standardized to have mean zero.
The assumption also holds if one of the covariates is a two-way interaction
term of two covariates that are tested in the two-stage strategy and are
included in the model as in the situation described in the main part of the
paper.

## Appendix C: FWER control

Let $T_{j}^{S_{1}},j=1,\ldots,p$ be the test-statistics for the $p$ covariates
in the first stage, and let $\Gamma_{j}^{S_{1}},j=1,\ldots,p$ be the
corresponding rejection regions. Furthermore, let $T_{jk}^{S_{2}}$ be the
test-statistic for testing interaction between the covariates $j$ and $k$ in
the second stage, and let $\Gamma_{jk}^{S_{2}}$ be the corresponding rejection
regions. In the second stage a null hypothesis is rejected if the p-value is
below the threshold $\alpha$ divided by the number of tests that are performed
in the second stage. This number is stochastic, and thus the rejection region
$\Gamma_{jk}^{S_{2}}$ is stochastic as well.

######  Lemma .5 (FWER control).

Suppose that, under the full null hypothesis, $(T_{j}^{S_{1}},T_{k}^{S_{1}})$
and $T_{jk}^{S_{2}}$ are asymptotically independent, and that, moreover,
$T_{j}^{S_{1}},j=1,\ldots,p$ are asymptotically independent. Then, under the
full null hypothesis, the two-step procedure preserves the FWER at the
significance level $\alpha$ in weak sense:

$\displaystyle\mbox{P}\Big{(}\bigcup_{j,k=1\mathrel{\mathop{\mathchar
58\relax}}j<k}^{p,p}\big{\\{}\big{(}T_{j}^{S_{1}}\in\Gamma_{j}^{S_{1}}\big{)}\cap\big{(}T_{k}^{S_{1}}\in\Gamma_{k}^{S_{1}}\big{)}\cap\big{(}T_{jk}^{S_{2}}\in\Gamma_{jk}^{S_{2}}\big{)}\big{\\}}\Big{)}\leq\alpha.$

######  Proof .6.

Let $R_{1}$ be defined as the stochastic number of tests that are rejected in
the first stage. Given that $R_{1}=r_{1}$ covariates were past on to the
second stage, the set $\Gamma_{jk}^{S_{2}}(r_{1})$ is defined as the non-
stochastic rejection region for test-statistic $T^{S_{2}}_{jk}$.

$\displaystyle\mbox{P}\Big{(}\bigcup_{j,k=1\mathrel{\mathop{\mathchar
58\relax}}j<k}^{p,p}\big{\\{}\big{(}T_{j}^{S_{1}}\in\Gamma_{j}^{S_{1}}\big{)}\cap\big{(}T_{k}^{S_{1}}\in\Gamma_{k}^{S_{1}}\big{)}\cap\big{(}T_{jk}^{S_{2}}\in\Gamma_{jk}^{S_{2}}\big{)}\big{\\}}\Big{)}$
$\displaystyle=\sum_{r_{1}=2}^{p}\mbox{P}(R_{1}=r_{1})\mbox{P}\Big{(}\bigcup_{j,k=1\mathrel{\mathop{\mathchar
58\relax}}j<k}^{p,p}\big{\\{}\big{(}T_{j}^{S_{1}}\in\Gamma_{j}^{S_{1}}\big{)}\cap\big{(}T_{k}^{S_{1}}\in\Gamma_{k}^{S_{1}}\big{)}\cap\big{(}T_{jk}^{S_{2}}\in\Gamma_{jk}^{S_{2}}\big{)}\big{\\}}\big{|}R_{1}=r_{1}\Big{)}$
$\displaystyle\leq\sum_{r_{1}=2}^{p}\mbox{P}(R_{1}=r_{1})\sum_{j,k=1\mathrel{\mathop{\mathchar
58\relax}}j<k}^{p,p}\mbox{P}\big{\\{}\big{(}T_{j}^{S_{1}}\in\Gamma_{j}^{S_{1}}\big{)}\cap\big{(}T_{k}^{S_{1}}\in\Gamma_{k}^{S_{1}}\big{)}\cap\big{(}T_{jk}^{S_{2}}\in\Gamma_{kj}^{S_{2}}\big{)}\big{|}R_{1}=r_{1}\big{\\}}$
$\displaystyle=\sum_{r_{1}=2}^{p}\mbox{P}(R_{1}=r_{1})\sum_{j,k=1\mathrel{\mathop{\mathchar
58\relax}}j<k}^{p,p}\mbox{P}\big{\\{}\big{(}T_{j}^{S_{1}}\in\Gamma_{j}^{S_{1}}\big{)}\cap\big{(}T_{k}^{S_{1}}\in\Gamma_{k}^{S_{1}}\big{)}\big{|}R_{1}=r_{1}\big{\\}}\mbox{P}\big{\\{}T_{jk}^{S_{2}}\in\Gamma_{jk}^{S_{2}}\big{|}\big{(}T_{j}^{S_{1}}\in\Gamma_{j}^{S_{1}}\big{)}\cap\big{(}T_{k}^{S_{1}}\in\Gamma_{k}^{S_{1}}\big{)}\cap\big{(}R_{1}=r_{1}\big{)}\big{\\}}$
$\displaystyle=\sum_{r_{1}=2}^{p}\mbox{P}(R_{1}=r_{1})\sum_{j,k=1\mathrel{\mathop{\mathchar
58\relax}}j<k}^{p,p}\mbox{P}\big{\\{}\big{(}T_{j}^{S_{1}}\in\Gamma_{j}^{S_{1}}\big{)}\cap\big{(}T_{k}^{S_{1}}\in\Gamma_{k}^{S_{1}}\big{)}\big{|}R_{1}=r_{1}\big{\\}}\mbox{P}\big{\\{}T_{jk}^{S_{2}}\in\Gamma_{jk}^{S_{2}}(r_{1})\big{|}\big{(}T_{j}^{S_{1}}\in\Gamma_{j}^{S_{1}}\big{)}\cap\big{(}T_{k}^{S_{1}}\in\Gamma_{k}^{S_{1}}\big{)}\cap\big{(}R_{1}=r_{1}\big{)}\big{\\}}$
$\displaystyle=\sum_{r_{1}=2}^{p}\mbox{P}(R_{1}=r_{1})\sum_{j,k=1\mathrel{\mathop{\mathchar
58\relax}}j<k}^{p,p}\mbox{P}\big{\\{}\big{(}T_{j}^{S_{1}}\in\Gamma_{j}^{S_{1}}\big{)}\cap\big{(}T_{k}^{S_{1}}\in\Gamma_{k}^{S_{1}}\big{)}\big{|}R_{1}=r_{1}\big{\\}}\mbox{P}\big{\\{}T_{jk}^{S_{2}}\in\Gamma_{jk}^{S_{2}}(r_{1})\big{\\}}$
$\displaystyle=\sum_{r_{1}=2}^{p}\mbox{P}(R_{1}=r_{1})\sum_{j,k=1\mathrel{\mathop{\mathchar
58\relax}}j<k}^{p,p}\mbox{P}\big{\\{}\big{(}T_{j}^{S_{1}}\in\Gamma_{j}^{S_{1}}\big{)}\cap\big{(}T_{k}^{S_{1}}\in\Gamma_{k}^{S_{1}}\big{)}\big{|}R_{1}=r_{1}\big{\\}}\frac{2\alpha}{r_{1}(r_{1}-1)}$
$\displaystyle=\sum_{r_{1}=2}^{p}\mbox{P}(R_{1}=r_{1})\sum_{j,k=1\mathrel{\mathop{\mathchar
58\relax}}j<k}^{p,p}\frac{r_{1}(r_{1}-1)}{p(p-1)}~{}~{}\frac{2\alpha}{r_{1}(r_{1}-1)}$
$\displaystyle=\;\alpha\sum_{r_{1}=2}^{p}\mbox{P}(R_{1}=r_{1})\;\leq\;\alpha.$

The fifth (in)equality follows from the asymptotic pairwise independence of
the test-statistics in the first and second stage.

In the theorem it is assumed that, under the full null hypothesis, the vector
$(T_{j}^{S_{1}},T_{k}^{S_{1}})$ and $T_{jk}^{S_{2}}$ are asymptotically
independent. This assumption follows from the asymptotic pairwise independence
proven in Result 2.2, and the fact that the test-statistics $T_{j}^{S_{1}}$
and $T_{k}^{S_{1}}$ are asymptotically Gaussian under the null hypothesis.

What is left, is to study the setting in which the test-statistics
$T_{j}^{S_{1}}$ and $T_{k}^{S_{1}}$ for any $j,k,j\neq k$ are asymptotically
independent. Define the models

$\displaystyle M_{1}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{X},{\beta}_{1},{\nu}_{1})\big{)}$
$\displaystyle=\beta_{1,1}x_{1}+\beta_{1,2}x_{2}+\beta_{1,3}x_{j}$
$\displaystyle M_{2}\mathrel{\mathop{\mathchar 58\relax}}\qquad
h\big{(}\mathbb{E}(Y|{X},{\beta}_{2},{\nu}_{2})\big{)}$
$\displaystyle=\beta_{2,1}x_{1}+\beta_{2,2}x_{2}+~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\beta_{2,4}x_{k}$

where $x_{1}=1$ for the intercept and $x_{2}$ the "fixed" covariate as before.
The situation is exactly the same as for proving complete independence. It
follows that for the Poisson, logistic and the linear regression model
asymptotic independence holds if the covariates that are included in the two
stage procedure (so not necessarily those that are permanently included in the
model: $x_{1}$ and $x_{2}$) are mutually independent and independent of the
“fixed” covariates, and the covariates are standardized to have zero mean,
independence of the test statistics $T_{j}^{S_{1}}$ and $T_{k}^{S_{1}}$ based
on the models $M_{1}$ and $M_{2}$ follows.

## Appendix D: Pairwise independence in the Cox PH model

In this Appendix the proof for pairwise independence of the regression
parameters in the Cox model is given, The proof is in line with the proof of
Theorem 2.1 in Subsection Appendix A, with the partial log likelihood function
in stead of the full log likelihood. However, since some nice properties of
the family of GLMs with canonical link functions do not hold for the Cox
model, the theorem only holds under the assumption of no main and interaction
effects.

Theorem 3.1 (Pairwise independence in the Cox PH model).  Let
$(T_{1}=\tilde{T}_{1}\wedge
C_{1},\Delta_{1},{X}_{1}),\ldots,(T_{n}=\tilde{T}_{n}\wedge
C_{1},\Delta_{n},{X}_{n})$ be a random sample. Suppose that
$\tilde{T}_{i}|X_{i}$ is a random observation from the Cox PH model with
hazard function as in (17), then $\sqrt{n}(\hat{\beta}_{1}-{\beta}_{1}^{0})$
and $\sqrt{n}\hat{\beta}_{2}^{\perp}$ are asymptotically independent.

######  Proof .7.

The partial log likelihood for the submodel $M_{1}$ in (15) is given by

$\displaystyle\frac{1}{n}PL_{1}({\beta}_{1})=\frac{1}{n}PL(E_{1}\beta_{1})=\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{\infty}\log\Bigg(\frac{\exp({X}_{i}^{t}E_{1}{\beta}_{1})}{\sum_{j=1}^{n}Y_{j}(t)\exp({X}_{j}^{t}E_{1}{\beta}_{1})}\Bigg{missing}){\rm
d}N_{i}(t),$

For ease of notation, we define $S^{(0)}$ and its first and second derivative
with respect to ${\beta}$, as

$\displaystyle S^{(0)}({\beta},t)$ $\displaystyle\mathrel{\mathop{\mathchar
58\relax}}=\;\frac{1}{n}\sum_{i=1}^{n}Y_{i}(t)\exp({X}_{i}^{t}{\beta})$
$\displaystyle S^{(1)}({\beta},t)$ $\displaystyle\mathrel{\mathop{\mathchar
58\relax}}=\;\frac{1}{n}\sum_{i=1}^{n}Y_{i}(t){X}_{i}\exp({X}_{i}^{t}\beta)$
$\displaystyle S^{(2)}({\beta},t)$ $\displaystyle\mathrel{\mathop{\mathchar
58\relax}}=\;\frac{1}{n}\sum_{i=1}^{n}Y_{i}(t){X}_{i}\exp({X}_{i}^{t}{\beta}){X}_{i}^{t}$

Then,

$\displaystyle\frac{1}{n}PL_{1}({\beta}_{1})=\frac{1}{n}PL(E_{1}{\beta}_{1})=\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{\infty}\log\Bigg(\frac{n^{-1}\exp({X}_{i}^{t}E_{1}{\beta}_{1})}{S^{(0)}(E_{1}{\beta}_{1},t)}\Bigg{missing}){\rm
d}N_{i}(t).$

Define $\hat{h}_{1}=\hat{\beta}_{1}-\beta_{1}^{0}$. By a Taylor expansion of
$PL_{1}$ with respect to $\hat{h}_{1}$ around zero

$\displaystyle\frac{1}{n}\widetilde{PL}_{1}(\hat{h}_{1})$
$\displaystyle\mathrel{\mathop{\mathchar
58\relax}}=\frac{1}{n}PL_{1}\big{(}E_{1}({\beta}_{1}^{0}+\hat{h}_{1})\big{)}$
$\displaystyle\propto\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{\infty}{X}_{i}^{t}E_{1}({\beta}_{1}^{0}+\hat{h}_{1})-\log\big(S^{(0)}(E_{1}({\beta}_{1}^{0}+\hat{h}_{1}),t)\big{missing})\;{\rm
d}N_{i}(t)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{\infty}{X}_{i}^{t}E_{1}({\beta}_{1}^{0}+\hat{h}_{1})-\log\big(S^{(0)}(E_{1}{\beta}_{1}^{0},t)\big{missing})-\hat{h}_{1}^{t}E_{1}^{t}~{}\frac{S^{(1)}(E_{1}{\beta}_{1}^{0},t)}{S^{(0)}(E_{1}{\beta}_{1}^{0},t)}$
$\displaystyle\qquad-\frac{1}{2}\hat{h}_{1}^{t}E_{1}^{t}\Big{(}\frac{S^{(2)}(E_{1}\beta_{1}^{0},t)}{S^{(0)}(E_{1}\beta_{1}^{0},t)}\;-\;\frac{S^{(1)}(E_{1}\beta_{1}^{0},t)^{t}S^{(1)}(E_{1}\beta_{1}^{0},t)}{S^{(0)}(E_{1}\beta_{1}^{0},t)^{2}}\Big{)}E_{1}\hat{h}_{1}\;{\rm
d}N_{i}(t)+O_{p}(\|\hat{h}_{1}\|^{3})$
$\displaystyle=K_{n}+\hat{h}_{1}^{t}E_{1}^{t}A_{n}+\frac{1}{2}\hat{h}_{1}^{t}E_{1}^{t}B_{n}E_{1}\hat{h}_{1}+O_{p}(\|\hat{h}_{1}\|^{3})$

with

$\displaystyle K_{n}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{\infty}{X}_{i}^{t}E_{1}{\beta}_{1}^{0}-\log\big(S^{(0)}(E_{1}\beta_{1}^{0},t)\big{missing}){\rm
d}N_{i}(t)$ $\displaystyle A_{n}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{\infty}{X}_{i}-\frac{S^{(1)}(E_{1}{\beta}_{1}^{0},t)}{S^{(0)}(E_{1}{\beta}_{1}^{0},t)}{\rm
d}N_{i}(t)$ $\displaystyle B_{n}$
$\displaystyle=-\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{\infty}\Big{(}\frac{S^{(2)}(E_{1}{\beta}_{1}^{0},t)}{S^{(0)}(E_{1}{\beta}_{1}^{0},t)}\;-\;\frac{S^{(1)}(E_{1}{\beta}_{1}^{0},t)^{t}S^{(1)}(E_{1}{\beta}_{1}^{0},t)}{S^{(0)}(E_{1}{\beta}_{1}^{0},t)^{2}}\Big{)}{\rm
d}N_{i}(t)$

The term $K_{n}$ does not depend on $\hat{h}_{1}$. Furthermore,
$\hat{h}_{1}^{t}E_{1}^{t}A_{n}=n^{-1/2}\hat{h}_{1}^{t}E_{1}^{t}n^{1/2}A_{n}$
where $n^{1/2}A_{n}$ equals $n^{-1/2}$ times the partial likelihood score
function in $\beta_{1}^{0}$, and, by the martingale central limit theorem (see
e.g., Andersen et al (1993)), is asymptotically normal with mean zero and
variance matrix denoted as $I$. That means that

$\displaystyle\hat{h}_{1}^{t}E_{1}^{t}A_{n}=\frac{1}{\sqrt{n}}\hat{h}_{1}^{t}E_{1}^{t}W_{1}+O_{p}(1/\sqrt{n})$

with $W_{1}$ normally distributed with mean zero and covariance matrix $I$.
The last term $B_{n}$ converges in probability to the covariance matrix $I$
(see e.g., Andersen et al (1993)). That means that

$\displaystyle\hat{h}_{1}^{t}E_{1}^{t}B_{n}E_{1}\hat{h}_{1}=\hat{h}_{1}^{t}E_{1}^{t}IE_{1}\hat{h}_{1}+o_{p}(\|\hat{h}_{1}\|^{2})$

The same can be done for submodel $M_{2}$, yielding the variable $W_{2}$.
Under the null hypothesis of no main and interaction effects, the covariance
of $W_{1}$ and $W_{2}$ is equal to
$\mathbb{E}W_{1}W_{2}^{t}=\mathbb{E}W_{1}W_{1}^{t}=I$.

The remainder of the proof follows the proof of Theorem 2.1 from Equation (22)
onwards.

## Appendix E: Power-plots for the Poisson regression model

In this appendix we show the power-plots for the Poisson regression model for
uncorrelated and correlated markers. The simulation study was performed as is
described in subsection 4.1 and 4.3.

Power as a function of the interaction effect for different first stage
threshold (FST) for the Poisson regression model. First row: independent
markers; Second row: correlated markers. First column: no main effects; Second
column: both main effects equal to 0.5; Third column: both main effects equal
to -0.5. (Note the different scales on the x-axes.)