Biometrics74, 1014–1022DOI: 10.1111/biom.12838

September 2018

Regression Analysis for Secondary Response Variable

in a Case-Cohort Study

Yinghao Pan

,

1

Jianwen Cai,

1

Sangmi Kim,

2

and Haibo Zhou

1,

*

1

Department of Biostatistics, University of North Carolina at Chapel Hill, Chapel Hill,

North Carolina 27599, U.S.A.

2

Medical College of Georgia, GRU Cancer Center, Augusta University, Augusta, Georgia 30912, U.S.A.

∗

email:[email protected]

Summary. Case-cohort study design has been widely used for its cost-effectiveness. In any real study, there are always other

important outcomes of interest beside the failure time that the original case-cohort study is based on. How to utilize the

available case-cohort data to study the relationship of a secondary outcome with the primary exposure obtained through the

case-cohort study is not well studied. In this article, we propose a non-parametric estimated likelihood approach for analyzing

a secondary outcome in a case-cohort study. The estimation is based on maximizing a semiparametric likelihood function

that is built jointly on both time-to-failure outcome and the secondary outcome. The proposed estimator is shown to be

consistent, efficient, and asymptotically normal. Finite sample performance is evaluated via simulation studies. Data from the

Sister Study is analyzed to illustrate our method.

Key words:Case-cohort design; Estimated likelihood; Secondary outcome; Semiparametric; Validation sample.

1. Introduction

The case-cohort study design (Prentice, 1986) is a cost-

effective sampling strategy that is used to study the

association between an often expensive exposure variable and

a time-to-event outcome. One example of a case-cohort study

is the Sister Study conducted by NIEHS (Kim et al., 2011,

2013), where the time-to-event outcome of interest is time to

breast cancer and the expensive exposure measure is a major

prostaglandin E2 metabolite (PGE-M). The Sister Study tar-

gets U.S. women who have a sister with breast cancer, but

with no breast cancer themselves at the enrollment. The case-

cohort study design can be viewed as a two-phased study. In

the first phase, information on the time-to-event outcome and

some relatively easy-to-obtain covariates are measured on all

cohort members. In the second phase, the exposure variable

of interest is measured in a random sample of the full cohort,

plus those subjects who experienced the event of interest. The

simple random sample here is referred as subcohort and those

failures are referred as the cases. Many methods have been

proposed to estimate the hazard ratio parameters with case-

cohort data, for example, Prentice (1986), Self and Prentice

(1988), Chen and Lo (1999), Borgan et al. (2000), Kulich and

Lin (2004), Kang and Cai (2009), Kim et al. (2013). For binary

response outcome-dependent-sampling data, related work was

done by Wang and Zhou (2006, 2010).

With tremendous cost and efforts involved in collecting

the exposure variable in case-cohort study, interest rises to

use the collected case-cohort data to study the relationship

between other important responses to the exposure. Relative

to the failure time used to design the case-cohort study, these

important responses are referred as the secondary responses.

For instance, in our Sister Study, investigators are interested

in studying the relationship between BMI and PGE-M since

recent research has indicated that there might be a positive

association between obesity and prostaglandin E2 (Morris

et al., 2011; Subbaramaiah et al., 2011). How to analyze the

secondary outcome (BMI in this case) in a case-cohort study

efficiently and correctly is not a straightforward exercise, as

the original study data are obtained in an outcome-dependent

way (depend on the primary outcome, i.e., the time-to-event

variable). This is an issue that puzzled investigators who try

to take advantage of the exposure variable measured in the

original case-cohort study, yet not sure how to handle the

biased sampling nature of the data based on the primary time-

to-event outcome. Naively treating the case-cohort sample as

a simple random sample could cause biased inference, as later

shown in our simulation studies. An inefficient approach is to

take the subcohort portion of the biased sampling data, and

ignore those of cases. Clearly, this approach is discarding a big

portion of the data. A significant amount of work was done on

secondary analysis of case-control or nested case-control data.

This includes the likelihood-based approach (Lee et al., 1997;

Jiang et al., 2006; Saarela et al., 2008; Lin and Zeng, 2009;

Salim et al., 2014); inverse probability weighting (Richardson

et al., 2007; Monsees et al., 2009); and estimating equation

(Wei et al., 2013; Ma and Carroll, 2016). Saarela et al. (2012)

proposed a conditional likelihood method for secondary anal-

ysis under general two-phase cohort sampling designs, which

includes case-cohort design. However, there has been a lack of

research on the secondary response regression analysis in the

case-cohort study in general.

Our research is motivated by the need in the Sister Study,

where we aim to establish the relationship between PGE-M

and BMI. In this article, we propose an estimated likelihood

1014©

2017, The International Biometric Society

Regression Analysis for Secondary Response Variable in a Case-Cohort Study1015

method for linear regression analysis of a continuous sec-

ondary response variable using case-cohort data. We jointly

model the time-to-event outcome (time to breast cancer)

and the continuous secondary response (BMI). The likeli-

hood function involves the conditional distribution of the

expensive exposure given other inexpensive covariates. We

estimate it in a non-parametric fashion. We compare our

proposed estimated likelihood estimator to two inverse prob-

ability weighting (IPW) type estimators we developed, and

show that the estimated likelihood method has greater sta-

tistical efficiency. The advantage of our proposed method

is that it is efficient, and yet require no strong parametric

assumptions. The performance of our estimator is explored

under a variety of conditions where complications could

arrive.

The organization of the article is as follows. In Section 2, we

present some notations, data structure, and model for case-

cohort design. In Section 3, we outline the estimation algo-

rithm for our proposed estimated likelihood estimator and

establish its asymptotic properties. We further develop two

new IPW type estimators in Section 4. In Section 5, we

investigate the finite sample performance of our proposed esti-

mators via simulation studies. In Section 6, we apply our

method to Sister Study data. Final remarks are given in

Section 7.

2. Data Structure and Model

We consider efficient inference of a continuous secondary

response, denoted byY

2

, with respect to an expensive expo-

sure, denoted byX, in a case-cohort design. To fix notation,

let

̃

Tdenote the primary event time of interest andCthe cen-

soring time. LetT=min(

̃

T,C), and=I(

̃

T≤C), whereTis

the observation time, andis the event indicator. Through-

out the article, we refer to individuals who have the event as

cases (=1) and censored individuals as non-cases (=0).

Furthermore, let (X,Z) denote the vector of covariates with

Xbeing the expensive scalar covariate obtained only for the

subcohort and the cases, andZbeing the other first-phase

covariates.Zcan be either discrete or continuous variables.

We assume that

̃

TandCare conditionally independent given

(Y

2

,X,Z), and the censoring timeCdoes not depend onX

but can depend onY

2

andZ.

We assume that the underlying data{(T

i

,

i

,Y

2i

,X

i

,Z

i

),i=

1,...,N}are independent and identically distributed random

vectors, whereNdenotes the size of the full cohort. Case-

cohort studies can be considered as two-phase studies: in the

first phase, information on observation time, event indicator,

secondary response, and inexpensive covariates are gathered

for each member of the full cohort. That is, we observe

{(T

i

,

i

,Y

2i

,Z

i

),i=1,...,N}. In the second phase, covariate

Xis measured for subjects in the subcohort, and those who

experienced the event of interest (

i

=1). LetV

0

be the index

set of the simple random sample taken from the baseline

cohort, andV

1

be the index set of the remaining cases in

the full cohort. Then we observe{X

i

,i∈V

0

∪V

1

}in the sec-

ond phase. Letn

SRS

be the size of the subcohort,n

V

1

be the

size of the supplemental cases. Here,n

SRS

is a pre-specified

number andn

V

1

is a random variable.

Letn

V

=n

SRS

+n

V

1

be the number of individuals for which

we observedX, and letn

̄

V

=N−n

V

be the number of

individuals for whomXis not observed. Borrowing the ter-

minology from measurement error literature, we will refer

to then

V

complete observations as the validation sample,

andn

̄

V

incomplete observations as non-validation sample. Let

V=V

0

∪V

1

be the index set of the validation sample, and

̄

Vbe the index set of the non-validation sample. In addi-

tion, letR

i

indicates whether theith subject is selected into

the validation sample, thenV={i:R

i

=1}and

̄

V={i:R

i

=

0}. The data structure for the problem we consider is the

following:

First phase:{T

i

,

i

,Y

2i

,Z

i

},i=1,...,N;

Second phase:{subcohort sample}{X

i

},i∈V

0

;

{supplemental cases}{X

i

|

i

=1,i /∈V

0

},i∈V

1

.

(1)

The main interest in this article is to model the association

betweenXand the secondary responseY

2

, adjusting forZin

the population from the following linear model:

Y

2

=β

0

+β

1

X+β

2

Z+,(2)

where∼N(0,σ

2

), and’s are independent. We note that

model (2) above, not model (3) below, is the main model

for this article. Our goal is to develop an efficient inference

procedure onβ=(β

0

,β

1

,β

2

).

Without loss of generality, we assume the following Cox

model (Cox, 1972) for the primary response variable:

λ(t|X, Y

2

,Z)=λ

0

(t) exp(γ

1

X+γ

2

Y

2

+γ

3

Z),(3)

whereλ(t|X, Y

2

,Z) is the conditional hazard function of

̃

T

given (X, Y

2

,Z), andλ

0

(t) is the baseline hazard function.

Let

0

(t)=

∫

t

0

λ

0

(u)dudenote the baseline cumulative hazard

function.

3. Estimated Likelihood Approach

Letf(·|·) denote the conditional density function, and

G(·|·) denote the conditional distribution function. LetW

be the informative components ofZaboutX, in the sense

thatG

X|Z

(x|z)=G

X|W

(x|w). If we are to jointly model

(T,  , Y

2

), then the full-information likelihood for data in (1)

is proportional to

L

N

(ξ,

0

(t))=

∏

i∈V

f(T

i

,

i

,Y

2i

|X

i

,Z

i

)

×

∏

j∈

̄

V

{

∫

x

f(T

j

,

j

,Y

2j

|x,Z

j

)dG

X|W

(x|W

j

)

}

,

(4)

1016Biometrics, September2018

whereξ=(β,γ,σ

2

). The log-likelihood of (4) takes the form:

l

N

(ξ,

0

(t))=

∑

i∈V

{

logf

γ

(T

i

,

i

|X

i

,Y

2i

,Z

i

)

+logf

β,σ

2

(Y

2i

|X

i

,Z

i

)

}

+

∑

j∈

̄

V

log

{

∫

x

f

γ

(T

j

,

j

|x, Y

2j

,Z

j

)

×f

β,σ

2

(Y

2j

|x,Z

j

)dG

X|W

(x|W

j

)

}

.(5)

The parameters contained in the above equation involve

(β,γ,σ

2

), though our main focus is on the inference ofβ

whileγ,σ

2

are nuisance parameters. Later, we will show that

our proposed estimated likelihood estimator is more efficient

than IPW type estimators because we took a joint likelihood

approach and incorporated the failure outcome information

(T, ) in the full cohort. Note that the conditional distribution

ofXgivenW,G

X|W

is involved in (4) and (5).

Thus, we propose to work with the following estimated

likelihood:

ˆ

l

N

(ξ,

0

(t))=

∑

i∈V

{

logf

γ

(T

i

,

i

|X

i

,Y

2i

,Z

i

)

+logf

β,σ

2

(Y

2i

|X

i

,Z

i

)

}

+

∑

j∈

̄

V

log

{

∫

x

f

γ

(T

j

,

j

|x, Y

2j

,Z

j

)

×f

β,σ

2

(Y

2j

|x,Z

j

)d

ˆ

G

X|W

(x|W

j

)

}

,(6)

where we non-parametrically estimate the conditional distri-

butionG

X|W

. For discreteW, let

ˆ

G

X|W

(x|w)=

∑

i∈V

0

I(X

i

≤x,W

i

=w)

/

∑

i∈V

0

I(W

i

=w).

For continuousW, we use the kernel method to estimate

G

X|W

.

ˆ

G

X|W

(x|w)=

∑

i∈V

0

I(X

i

≤x)K

H

(W

i

−w)

/

∑

i∈V

0

K

H

(W

i

−w),

whereK

H

(·)=|H|

−1/2

K(H

−1/2

·) is a multivariate kernel with

a bandwidth matrixH. The optimal bandwidthHcan be

selected using cross-validation, or using an ad hoc value. For

instance, whenWis univariate, a simple bandwidth selec-

tion isH=2

̂σ

W

·n

−1/3

SRS

, wherê

σ

W

is the estimated standard

deviation ofW.

In the estimated likelihood expression (6),f

γ

(T

i

,

i

|

X

i

,Y

2i

,Z

i

) andf

β,σ

2

(Y

2i

|X

i

,Z

i

) can be written in explicit

forms based on linear model (2) and Cox model (3). Under

the assumption that

̃

TandCare conditionally independent,

f

γ

(T

i

,

i

|X

i

,Y

2i

,Z

i

) is equal to

{

λ

0

(T

i

) exp(γ

1

X

i

+γ

2

Y

2i

+γ

3

Z

i

)

}



i

e

−

0

(T

i

) exp(γ

1

X

i

+γ

2

Y

2i

+γ

3

Z

i

)

λ

1−

i

C

(T

i

|X

i

,Y

2i

,Z

i

)e

−

C

(T

i

|X

i

,Y

2i

,Z

i

)

,

whereλ

C

(T|X, Y

2

,Z) is the hazard function ofCgiven

(X, Y

2

,Z), and

C

is the cumulative hazard. Assuming that

Cis independent ofX, logf

γ

(T

i

,

i

|X

i

,Y

2i

,Z

i

) is the sum of



i

logλ

0

(T

i

)+

i

(γ

1

X

i

+γ

2

Y

2i

+γ

3

Z

i

)−

0

(T

i

) exp(γ

1

X

i

+γ

2

Y

2i

+γ

3

Z

i

)

and (1−

i

) logλ

C

(T

i

|Y

2i

,Z

i

)−

C

(T

i

|Y

2i

,Z

i

). The latter term

does not involvexor (ξ,

0

(t)), and thus irrelevant to estimat-

ing the parameters of interest. Also,

logf

β,σ

2

(Y

2i

|X

i

,Z

i

)=−

1

2

log(2

πσ

2

)−

1

2σ

2

(Y

2i

−β

0

−β

1

X

i

−β

2

Z

i

)

2

.

The estimated log-likelihood

ˆ

l

N

(ξ,

0

(t)) has a semi-

parametric form, withξbeing the parametric part and

0

(t)

the non-parametric component. Because the baseline cumu-

lative hazard

0

(t) is unknown, we cannot directly maximize

ˆ

l

N

(ξ,

0

(t)) with respect toξ. We propose to use an inverse

probability weighted partial likelihood estimates from Chen

and Lo (1999) to obtain (

ˆ

γ

0

1

,

ˆ

γ

0

2

,

ˆ

γ

0

3

) as initial parameter esti-

mates for the Cox model (3). We then estimate the baseline

cumulative hazard using a weighted Breslow estimator,

ˆ



0

(t)=

∑

j∈V

I(T

j

≤t)

j

∑

l∈V

ˆ

π

−1

l

I(T

l

≥T

j

) exp(

ˆ

γ

0

1

X

l

+

ˆ

γ

0

2

Y

2l

+

ˆ

γ

0

3

Z

l

)

,

where

ˆ

π

l

=1 for cases, and

ˆ

π

l

=n

SRS

/Nfor non-cases. Fur-

thermore, we use the inverse probability weighted linear

regression on the validation sample to obtain (

ˆ

β

0

0

,

ˆ

β

0

1

,

ˆ

β

0

2

,

ˆ

σ

0

)

as initial parameter estimates for the linear regression model

(2), where cases receive unit weight, and non-cases have

weightN/n

SRS

.

Finally, we plug in the baseline cumulative hazard estimate

ˆ



0

(t) into the estimated log-likelihood, and maximize it with

respect toξ. Let

̃

l

N

(ξ)=

ˆ

l

N

(ξ,

ˆ



0

(t)),

then our proposed non-parametric estimated likelihood esti-

mator

ˆ

ξ

NPEL

is the solution to the following estimating

equation (7):

1

N

∂

̃

l

N

(ξ)

∂ξ

=0.(7)

Newton–Raphson algorithm can now be invoked, where

(

ˆ

γ

0

1

,

ˆ

γ

0

2

,

ˆ

γ

0

3

,

ˆ

β

0

0

,

ˆ

β

0

1

,

ˆ

β

0

2

,

ˆ

σ

0

) are the starting values.

Next, we present the asymptotic property of our proposed

estimator. The following results are provided for the case

whenWis discrete. Asymptotic results forWbeing con-

tinuous can be derived similarly. We assume that asN→∞,

n

V

/N→

̄

π>0, andn

SRS

/n

V

→ρ

SRS

>0. Therefore,

̄

πis the

proportion of validation sample in the overall population.ρ

SRS

represents the proportion of simple random sample in the val-

idation sample. In addition, we letE

k

denote the conditional

expectation given=k. That is,E

k

{h(T,  , Y

2

,X,Z)}=

E{h(T,  , Y

2

,X,Z)|=k}.Letξ

∗

be the true underlying

value ofξ. Under some general regularity conditions (see Web

Appendix), we have the following theorems.

Regression Analysis for Secondary Response Variable in a Case-Cohort Study1017

Theorem1.(i)

ˆ

ξ

NPEL

is a consistent estimator ofξ

∗

and

(ii)

ˆ

ξ

NPEL

has the following asymptotic distributional proper-

ties:

√

N(

ˆ

ξ

NPEL

−ξ

∗

)

D

−→N(0,

NPEL

(ξ

∗

)).

where

NPEL

(ξ

∗

)=I

−1

(ξ

∗

)+ρ

SRS

̄

πI

−1

(ξ

∗

)(ξ

∗

)I

−1

(ξ

∗

),(ξ)=

var

X,W

{

(1−

̄

π)Q

X,W

(ξ,

0

)

}

, and

I(ξ)=−ρ

SRS

̄

πE

{

∂

2

logf(T,  , Y

2

|X,Z;ξ)

∂ξ∂ξ

T

}

−(1−

̄

π)E

0

{

∂

2

logf(T,  , Y

2

|Z;ξ)

∂ξ∂ξ

T

}

−(1−ρ

SRS

)

̄

πE

1

{

∂

2

logf(T,  , Y

2

|X,Z;ξ)

∂ξ∂ξ

T

}

,(8)

Q

X,W

(ξ,

0

)=

∑

w

P(W=w|=0)

I(W=w)

P(W=w)ρ

SRS

̄

π

×E

T, ,Y

2

,Z|=0

{

M

X

(T,  , Y

2

,Z,W=w;ξ,

0

)

}

.

(9)

where

M

X

(T,  , Y

2

,Z,W;ξ,

0

)=

∂f

(T,  , Y

2

|X,Z;

0

)/∂ξ

∫

f(T,  , Y

2

|X,Z;

0

)dG

X|W

(X|W)

−

f

(T,  , Y

2

|X,Z;

0

)

∫

∂f(T,  , Y

2

|X,Z;

0

)/∂ξdG

X|W

(X|W)

|

∫

f(T,  , Y

2

|X,Z;

0

)dG

X|W

(X|W)|

2

.

(10)

Replacing the population quantities in Theorem 1 with

the sample quantities, a consistent estimator for asymptotic

covariance matrix

NPEL

(ξ

∗

) is provided by the following

theorem:

Theorem2.A consistent estimator for variance matrix



NPEL

(ξ

∗

)can be obtained as

̂



NPEL

(

ˆ

ξ

NPEL

)=

ˆ

I

−1

(

ˆ

ξ

NPEL

)+

n

SRS

N

×

ˆ

I

−1

(

ˆ

ξ

NPEL

)

̂

(

ˆ

ξ

NPEL

)

ˆ

I

−1

(

ˆ

ξ

NPEL

),(11)

where

ˆ

I(ξ)=−

1

N

∂

2

̃

l

N

(ξ)

∂ξ∂ξ

T

,and

̂

(ξ)=̂var

{X

i

,W

i

:i∈V

0

}

{

N−n

V

N

̂

Q

X

i

,W

i

(ξ)

}

,with

̂

Q

X

i

,W

i

(ξ)=

∑

j∈

̄

V

1

N−n

V

×

N·I(W

i

=W

j

)

∑

i∈V

0

I(W

i

=W

j

)

×

[

∂f(T

j

,

j

,Y

2j

|X

i

,Z

j

;

ˆ



0

)/∂ξ

ˆ

f(T

j

,

j

,Y

2j

|Z

j

;

ˆ



0

)

−

f(T

j

,

j

,Y

2j

|X

i

,Z

j

;

ˆ



0

)∂

ˆ

f(T

j

,

j

,Y

2j

|Z

j

;

ˆ



0

)/∂ξ

{

ˆ

f(T

j

,

j

,Y

2j

|Z

j

;

ˆ



0

)

}

2

]

,

where

ˆ

f(T

j

,

j

,Y

2j

|Z

j

;

ˆ



0

)=

∫

f(T

j

,

j

,Y

2j

|x,Z

j

;

ˆ



0

)d

ˆ

G

X|W

(x|W

j

).

Outline of the proofs of these two theorems are provided in

Web Appendix A.

4. Derivation of Two New IPW Type Estimators

In this section, we derive two inverse probability weight-

ing type estimators for the secondary outcome analysis. One

is the basic weighted estimator and the other one is the

augmented IPW estimator. Though the general ideas of prob-

ability weighted estimator is there (Horvitz and Thompson,

1952), we develop out these ideas explicitly as competing esti-

mators for the proposed

ˆ

ξ

NPEL

. Later in the simulation study,

we will compare our proposed estimator to these two IPW

type estimators.

The first IPW estimator

ˆ

β

IPW

, is a basic weighted least

squares estimator, where the weight is taken as the inverse of

the selection probability into the validation sample. That is,

we want to minimize

∑

i∈V

1

π

i

(Y

2i

−β

0

−β

1

X

i

−β

2

Z

i

)

2

=

N

∑

i=1

R

i

π

i

(Y

2i

−β

0

−β

1

X

i

−β

2

Z

i

)

2

,

whereπ

i

is the selection probability for theith subject. In the

case-cohort design,π

i

=1 for cases, as all cases are included

into the validation sample. For non-cases,π

i

=n

SRS

/N, the

probability of being selected into the SRS portion of the

validation sample. Let

̃

X

i

=(1,X

i

,Z

i

) be the vector of covari-

ates for theith subject, including the intercept. It is easy to

see that the IPW estimator

ˆ

β

IPW

satisfies the following score

equation:

N

∑

i=1

R

i

π

i

̃

X

T

i

(Y

2i

−

̃

X

i

β)=0.(12)

We notice that in (12), all cases will receive unit weights,

and all non-cases selected into the validation sample will have

weights larger than one. The inverse probability weighting

method will give an unbiased estimator of the parametersβ.

However it may also result in reduced efficiency, since the

cases would always receive smaller weights than non-cases

regardless of whether the case status is linked to the secondary

outcome/exposure relationship.

The second IPW estimator,

ˆ

β

AIPW

, is an augmented IPW

that incorporates the available information in the full cohort.

Similar to Robins et al. (1994), we have the following estimat-

ing equation:

N

∑

i=1

[

R

i

π

i

̃

X

T

i

(Y

2i

−

̃

X

i

β)+

(

1−

R

i

π

i

)

E

X

i

|Y

2i

,Z

i

{

̃

X

T

i

(Y

2i

−

̃

X

i

β)

}

]

=0,

(13)

where an augmented term is added. The second term in

(13) involvesE(X

i

|Y

2i

,Z

i

) andE(X

2

i

|Y

2i

,Z

i

). We use a lin-

ear regression to approximate these moments, that is,E(X

i

|

Y

2i

,Z

i

)=φ

0

+φ

1

Y

2i

+φ

2

Z

i

, Var(X

i

|Y

2i

,Z

i

)=τ

2

. Based on

1018Biometrics, September2018

SRS portion of the data, a linear regression can be fitted to

obtain parameter estimates for (φ

0

,φ

1

,φ

2

,τ

2

). After plugging

(

ˆ

φ

0

,

ˆ

φ

1

,

ˆ

φ

2

,

ˆ

τ

2

) into (13), Newton–Raphson algorithm can be

invoked to obtain

ˆ

β

AIPW

.

5. Simulation Studies

We conducted a series of simulation studies. First, we evalu-

ated the performance of our estimator when the underlying

assumption is met. Then, we further broaden out to evaluate

the proposed method under situations where models are mis-

specified: (i) IfW, the informative components ofZabout

X, is misspecified. (ii) The error termin linear model (2)

is not normally distributed. (iii) The proportional hazards

assumption does not hold.

Five competing estimators are compared: (i)

ˆ

β

N

, linear

regression of the secondary responseY

2

on (X,Z) based on

the validation sample (union of subcohort and cases). This

naive estimator ignores the biased sampling nature of the

data, and can be grossly misleading; (ii)

ˆ

β

SRS

, fitting a lin-

ear regression ofY

2

on (X,Z) based on the SRS portion of

the validation sample; (iii)

ˆ

β

IPW

denotes the inverse prob-

ability weighting method using only the observations from

the validation sample; (iv)

ˆ

β

AIPW

denotes the augmented

inverse probability weighting method we derive in Section

4 which incorporates the available information in the full

cohort; (v)

ˆ

β

NPEL

denotes the estimated likelihood method we

proposed.

In the first set of studies, we let the full cohort size

to beN=1000. The data were generated by the follow-

ing models:Y

2

=β

0

+β

1

X+β

2

Z+,whereX∼N(0,1),Z∼

Bernoulli(0.45), and∼N(0,1). We setβ

0

=1,β

2

=−0.5,

and allowβ

1

to take value 0 or 0.5. The event time

̃

T

follows an exponential distribution with hazard function:

λ(t|X, Y

2

,Z)=λ

0

(t) exp(γ

1

X+γ

2

Y

2

+γ

3

Z),withλ

0

(t)=1,

γ

1

=0.5,γ

3

=−0.5, andγ

2

takes value 0 or log(3). We assume

that the censoring time follows a uniform distribution on

interval (0,C). Whenγ

2

=0,C=0.52 or 0.15 corresponds

to 80% or 95% censoring, respectively. Whenγ

2

=log(3),

C=0.16 or 0.02 corresponds to 80% or 95% censoring,

respectively.

From the full cohort, we first select a simple random sample

of sizen

SRS

=200, then we add all cases for the case-cohort

sample. As we are mainly interested in studying the relation-

ship between the secondary outcomeY

2

and the expensive

exposureX, we focus on the estimation ofβ

1

. We assume that

the informative components ofWisZ. SinceZis discrete, we

use the empirical distribution function to estimate the condi-

tional distribution ofXgivenZ. We report the mean of the

parameter estimates, empirical standard deviation, mean of

the estimated standard deviation, and 95% confidence inter-

val coverage. Simulation results based on 1000 replications are

shown in Table 1.

We find that all estimators yield approximately unbiased

estimates except

ˆ

β

N

.

ˆ

β

N

performs well when the secondary

outcome is not correlated with the time-to-event outcome

(γ

2

=0). However, it has bias and poor coverage rate when

γ

2

=0. The bias is even larger when censoring percent

increases (rare event). Similar findings were reported by Lin

and Zeng (2009) in secondary analysis of case-control data.

Table 1

Simulation results with full cohort sizeN=1000, the SRS

sample size isn

SRS

=200.X∼N(0,1),Z∼Bernoulli(0.45).

Estimated quantities forβ

1

Trueβ

1

Trueγ

2

MethodsBiasSD

̂

SDCI

Censoring percentage = 80%

00

ˆ

β

N

0.001 0.050 0.052 0.950

ˆ

β

SRS

0.000 0.070 0.071 0.948

ˆ

β

IPW

0.002 0.063 0.063 0.950

ˆ

β

AIPW

0.002 0.064 0.063 0.944

ˆ

β

NPEL

0.002 0.051 0.052 0.945

log(3)

ˆ

β

N

0.024 0.054 0.054 0.925

ˆ

β

SRS

0.000 0.071 0.071 0.947

ˆ

β

IPW

0.001 0.062 0.061 0.947

ˆ

β

AIPW

0.002 0.063 0.062 0.944

ˆ

β

NPEL

0.001 0.053 0.052 0.945

0.50

ˆ

β

N

−0.003 0.051 0.052 0.959

ˆ

β

SRS

−0.005 0.071 0.071 0.953

ˆ

β

IPW

−0.005 0.063 0.063 0.939

ˆ

β

AIPW

−0.002 0.054 0.055 0.954

ˆ

β

NPEL

−0.005 0.046 0.048 0.955

log(3)

ˆ

β

N

0.025 0.052 0.051 0.922

ˆ

β

SRS

−0.001 0.071 0.071 0.945

ˆ

β

IPW

0.000 0.063 0.060 0.935

ˆ

β

AIPW

0.004 0.056 0.054 0.944

ˆ

β

NPEL

−0.001 0.047 0.047 0.950

Censoring percentage = 95%

0.5log(3)

ˆ

β

N

0.084 0.066 0.065 0.729

ˆ

β

SRS

−0.001 0.073 0.071 0.940

ˆ

β

IPW

0.002 0.069 0.066 0.928

ˆ

β

AIPW

0.005 0.062 0.059 0.933

ˆ

β

NPEL

0.000 0.058 0.058 0.943

Note: The models areY

2

=β

0

+β

1

X+β

2

Z+,λ(t|X, Y

2

,Z)=

λ

0

(t) exp(γ

1

X+γ

2

Y

2

+γ

3

Z).

From Table 1, we draw several other conclusions. (i) The

proposed estimated likelihood estimator

ˆ

β

NPEL

is the most

efficient in all settings. (ii) The estimated standard devia-

tion is very close to the empirical standard deviation (i.e.,

̂

SDis close to SD). (iii) The 95% confidence interval cover-

age is close to 0.95, which implies that the asymptotic normal

approximation works well in these situations. (iv) We also

notice that whenβ

1

is significantly different from 0 (i.e.,

β

1

=0.5),

ˆ

β

AIPW

outperforms

ˆ

β

IPW

by a certain margin. The

difference between these two estimators is negligible when

β

1

=0.

We conducted further simulation studies to assess the per-

formance of our proposed estimator under a wide range of

proportion of SRS sample in the full cohort, and under vari-

ous censoring percentages. Simulation results are provided in

Web Appendix B.

Regression Analysis for Secondary Response Variable in a Case-Cohort Study1019

5.1.Misspecification ofW, the Informative Components

ofZAboutX

To study the performance of

ˆ

β

NPEL

whenWis misspecified,

we generate the data by the following models:Y

2

=β

0

+β

1

X+

β

2

Z

1

+β

3

Z

2

+,whereZ

1

∼Bernoulli(0.45),Z

2

∼Bernoulli

(0.6),∼N(0,1). The expensive exposureXis simulated

as:X∼N(Z

1

,1). LetZ=(Z

1

,Z

2

), thenG(X|Z)=G(X|Z

1

),

that is, the trueWisZ

1

. We setβ

0

=1,β

2

=−0.5,β

3

=

2, and allowβ

1

to take value 0 or 0.5. The event

time

̃

Tfollows an exponential distribution with haz-

ard function:λ(t|X, Y

2

,Z)=λ

0

(t) exp(γ

1

X+γ

2

Y

2

+γ

3

Z

1

+

γ

4

Z

2

),withλ

0

(t)=1,γ

1

=0.5,γ

2

=log(1.2),γ

3

=−0.5 and

γ

4

=1. The censoring time follows a uniform distribution

on interval (0,0.14), which leads to 80% censoring. Same as

before,N=1000 andn

SRS

=200. We consider three scenar-

ios: (i) We estimate the conditional distributionG(X|Z)by

assuming thatWisZ

2

. This is the case whereWis misspec-

ified. (ii) We assume thatWisZ

1

. (iii) We assume thatW

is (Z

1

,Z

2

).

Table 2 summarizes the simulation results. Under null

hypothesis (β

1

=0),

ˆ

β

NPEL

has good performance even if we

misspecifyW. However, whenβ

1

is significantly different from

0 (i.e.,β

1

=0.5), misspecification ofW(usingZ

2

only) would

lead to substantial bias and poor coverage probability of 95%

CI. In the Discussion Section, we will provide practical guid-

ance on how to chooseW. The general idea is that we should

select the components ofZthat correlates toXthe most.

5.2.Parametric versus Non-Parametric Estimation ofX

GivenZ

We compare our proposed estimator

ˆ

β

NPEL

to a parametric

maximum likelihood estimator

ˆ

β

PL

, where we parametrically

specifyG

X|Z

(i.e., conditional normal) to a set of additional

parametersλ, and then directly maximize the log-likelihood

function in (5) with respect toξandλ. The simulation set up

is the same as Table 2, except forX. We consider the follow-

ing three scenarios: (i)X∼N(0,1). (ii)Xfollows a gamma

distribution with shape parameter 0.5 and rate parameter 1.

(iii)Xfollows an exponential distribution with rate parameter

0.5.

Table 2

Performance of our proposed estimator

ˆ

β

NPEL

when we

misspecifyW.Z=(Z

1

,Z

2

), whereZ

1

∼Bernoulli(0.45),

Z

2

∼Bernoulli(0.6).X∼N(Z

1

,1), which means thatZ

1

is

the informative components ofZaboutX. The trueWisZ

1

.

Estimated quantities forβ

1

Trueβ

1

WusedBiasSD

̂

SDCI

0Z

1

0.0010.0530.0520.952

Z

2

0.0010.0490.0490.936

(Z

1

,Z

2

)0.0010.0530.0520.954

0.5Z

1

−0.0020.0460.0470.958

Z

2

−0.0630.0440.0450.728

(Z

1

,Z

2

)−0.0070.0460.0480.964

Note: Models areY

2

=β

0

+β

1

X+β

2

Z

1

+β

3

Z

2

+,λ(t|X, Y

2

,Z)=

λ

0

(t) exp(γ

1

X+γ

2

Y

2

+γ

3

Z

1

+γ

4

Z

2

).

Table 3 displays the results.

ˆ

β

PL

has large bias and can be

grossly misleading whenG

X|Z

is misspecified. On the other

hand, our proposed estimator

ˆ

β

NPEL

works well in all settings.

Another disadvantage of parametric MLE is that it requires

numerical integration overX, and is hence computationally

intensive.

5.3.Misspecification of the Error Distribution for the

Secondary Outcome or Misspecification of the

Survival Model

We conducted additional simulations to examine the perfor-

mance of our estimator for the secondary outcome parameters

when the error distribution for the secondary outcome model

is misspecified or when the survival model is misspecified.

The simulation set up is the same as Table 2, except for the

error term and the model for time to event outcome. We con-

sider five scenarios: (i) the error term is a gamma distribution

with shape parameter 2, rate parameter 1, then normalized to

have mean 0 and variance 1. This error term is right-skewed.

(ii) The error term depends onZ

1

, which is a standard nor-

mal distribution multiplied by (1+Z

2

1

)

3/4

/2. (iii) The error

term depends onX, which is a standard normal distribution

multiplied by (1+X

2

)

3/4

/2. (iv) The event time

̃

Tcomes

from an Accelerated Failure Time (AFT) model: log

̃

T=

0.5X+Y

2

−0.5Z

1

+Z

2

+e, wheree∼N(0,1). (v) The event

time follows an Additive Hazards model:λ(t)=10+0.5X+

Y

2

−0.5Z

1

+Z

2

.

Simulation results are summarized in Table 4, where we use

Z

1

to estimateG(X|Z). We have the following conclusions:

(a) when the error term is skewed (not symmetric around 0),

ˆ

β

IPW

and

ˆ

β

AIPW

performs relatively well, while

ˆ

β

NPEL

has larger

bias. This comes from the fact that the estimated likelihood

estimator is more reliant on the normality assumption. On the

other hand, the IPW and AIPW estimator we developed do

not explicitly use the normality assumption of the error term.

(b) When the error term is heteroscedastic, but not depend

Table 3

Parametric versus non-parametric estimation ofXgivenZ

Estimated quantities forβ

1

MethodsBiasSD

̂

SDCI

Xis normal

ˆ

β

PL

0.0050.0490.0640.926

ˆ

β

NPEL

−0.0050.0500.0500.957

Xfollows a gamma

distribution

ˆ

β

PL

−0.0230.0550.0960.913

ˆ

β

NPEL

−0.0050.0590.0600.962

Xfollows an exponential

distribution

ˆ

β

PL

−0.0510.0170.0150.143

ˆ

β

NPEL

−0.0010.0170.0170.956

Note: Models areY

2

=β

0

+β

1

X+β

2

Z

1

+β

3

Z

2

+.λ(t|X, Y

2

,Z)=

λ

0

(t) exp(γ

1

X+γ

2

Y

2

+γ

3

Z

1

+γ

4

Z

2

).

ˆ

β

NPEL

,proposednon-

parametric estimated likelihood estimator.

ˆ

β

PL

, parametric

maximum likelihood estimator assuming thatXgivenZis normal.

1020Biometrics, September2018

Table 4

Simulation results when the error termis misspecified or the survival model is misspecified.Z=(Z

1

,Z

2

), where

Z

1

∼Bernoulli(0.45),Z

2

∼Bernoulli(0.6),X∼N(Z

1

,1).β

1

=0.5.

Estimated quantities forβ

1

MethodsBiasSD

̂

SDCI

Misspecification of

(i) Gamma error misspecified asN(0,σ

2

)

ˆ

β

IPW

−0.0010.0630.0610.936

ˆ

β

AIPW

0.0070.0650.0610.916

ˆ

β

NPEL

−0.0330.0510.0490.889

(ii) Heteroscedastic error (related toZ

1

) misspecified asN(0,σ

2

)

ˆ

β

IPW

0.0010.0440.0420.937

ˆ

β

AIPW

0.0040.0370.0360.942

ˆ

β

NPEL

−0.0040.0320.0320.956

(iii) Heteroscedastic error (related toX) misspecified asN(0,σ

2

)

ˆ

β

IPW

0.0010.0990.0960.940

ˆ

β

AIPW

0.0030.0960.0930.928

ˆ

β

NPEL

−0.0130.0740.0510.818

Misspecification of survival model

(iv) AFT misspecified as PH

ˆ

β

IPW

0.0050.0640.0630.940

ˆ

β

AIPW

0.0070.0570.0560.930

ˆ

β

NPEL

−0.0050.0490.0490.951

(v) Additive Hazards misspecified as PH

ˆ

β

IPW

0.0020.0670.0640.937

ˆ

β

AIPW

0.0080.0600.0560.912

ˆ

β

NPEL

0.0010.0490.0490.948

Note: Models for secondary outcome isY

2

=β

0

+β

1

X+β

2

Z

1

+β

3

Z

2

+. AFT, Accelerated Failure Time; PH, Proportional Hazards.

onX, the proposed

ˆ

β

NPEL

is unbiased and the most efficient

among all estimators. (c) When the error term depends onX,

ˆ

β

NPEL

is still unbiased, but the asymptotic covariance estima-

tor tends to underestimate the empirical covariance. (d) Even

though the true model for

̃

Tis AFT or Additive Hazards and

we fit Cox model, our proposed estimator

ˆ

β

NPEL

still works

very well in terms of estimating the parameters in model (2).

In practice, we could first fit a linear regression ofY

2

on

(X,Z) based on SRS portion of the data to check the normal-

ity and homogeneity assumption. If we find that the residual

depends onX, we need to employ some variable transforma-

tion techniques to obtain a more homogeneous error term.

After that, our estimated likelihood method can be used.

6. Analysis of Sister Study Data

We applied our method to data from the Sister Study (Kim

et al., 2011, 2013) to assess the relationship between BMI

and PGE-M. The Sister Study is a cohort study conducted

on U.S. women aged 35–74. To be eligible, women cannot have

breast cancer themselves at the start of the study but have

a sister with breast cancer. The participants come from all

50 states and Puerto Rico, and were recruited by health pro-

fessionals, internet, and a national campaign. The enrollment

time is from 2003 to 2009. Twenty-five thousand and eight

hundred subjects completed the baseline activities by June 1,

2007. We further restrict the study population to women who

are postmenopausal, age 50 or older, and not currently using

hormone replacement therapies, which yields a full cohort of

size 11,338.

Information on subject’s breast cancer status is followed-up

via annual and biennial questionnaires until September 2010.

Other information such as age, BMI, medication use are col-

lected for each member in the full cohort (N=11,338). In the

original study, the primary interest is to assess the effect of

PGE-M on the risk of breast cancer. Because the measure-

ment of PGE-M level from urine samples can be expensive,

a case-cohort design is implemented. A simple random sam-

ple of 300 participants were selected from the eligible full

cohort. Then 307 supplemental cases were added to create

a case-cohort sample. PGE-M levels were quantified for the

subcohort members and the supplemental cases.

Prior research indicate that obesity might be posi-

tively related to prostaglandin E2 (Morris et al., 2011;

Subbaramaiah et al., 2011). In this analysis, we use normal-

ized BMI as the secondary outcome and normalized PGE-M

as the exposure which is only available on the case-cohort

sample. Additional confounding variables include use of

Regression Analysis for Secondary Response Variable in a Case-Cohort Study1021

nonsteroidal anti-inflammatory drugs (NSAIDS), current

drinking status, current smoking status, indicator variables of

older than 60 and waist circumference larger than 35 inches.

We fit the following two models:

BMI=β

0

+β

1

PGEM+β

2

NSAIDS+β

3

DRINKING

+β

4

AGE+,

λ(t)=λ

0

(t) exp(γ

1

PGEM+γ

2

BMI+γ

3

NSAIDS+γ

4

AGE

+γ

5

WAIST+γ

6

SMOKE),

whereλ(t) is the conditional hazard of breast cancer given

PGE-M, BMI and other covariates. The estimated conditional

distribution of the expensive exposure PGE-M given other

covariates is based on use of NSAIDS, current drinking status,

and an indicator variable of older than 60.

The results are shown in Table 5. All three analyses confirm

that PGE-M level is significantly related to BMI. The pro-

posed estimator in general provides a more precise estimate

of the effects. Most of the variance estimates for the proposed

method are smaller than those from competing methods. All

three analyses confirm that there is a positive relationship

between BMI and use of NSAIDS, and that there is a negative

relationship between current drinking, age, and BMI.

7. Discussion

To study the hard-to-obtain exposure variable that is only

measured in a case-cohort study with some other important

responses is desired by all investigators. However, how to

properly account for the biased sampling nature of a case-

cohort study in a secondary outcome regression analysis, or

how to best utilize the available data to achieve a more effi-

cient inference is not a well thought out issue. Without loss

of generality, we assume a linear model for the secondary

response and our goal in this article is to obtain unbiased and

efficient estimates for the linear model parameters in model

(2). We propose a new and efficient estimator for the sec-

ondary outcome regression analysis in a case-cohort study.

Our estimation is based on jointly modeling the time-to-event

outcome and the secondary response. The joint likelihood

is not tractable as the conditional distribution of expensive

exposure given other inexpensive covariates as well as the

baseline cumulative hazard are not specified. We estimate this

conditional distribution and the baseline cumulative hazard

function non-parametrically. Our proposed method has the

advantage of making no new parametric assumptions, and

yet still able to improve efficiency relative to other competing

estimators.

In the simulation studies, we find that misspecification of

Wwould lead to certain degree of bias. We suggest the follow-

ing practical guidelines on how to chooseW. Due to “curse of

dimensionality,” theWwe use should not be more than three

components in general, unless the SRS sample size is very

large. When covariatesZhas low dimension (i.e.,≤3), one

can simply include all covariates intoW. When the dimen-

sion ofZis medium (i.e., 3< dim(Z)<10), we suggest to

calculate the correlation betweenXand each covariate in

Z, and then select three most correlated components. When

the dimension ofZis high (i.e.,≥10), we suggest to use

dimension-reduction techniques such as principal component

analysis to get a lower dimensional

̃

Z, then our proposed esti-

mated likelihood method can be carried out by using

̃

Zas the

regression covariates and the informative componentsW.

Our proposed estimated likelihood method requires spec-

ifying a model for the time-to-event outcome. However, in

the simulation studies, we found that our proposed estimator

works very well in estimating the parameters in the secondary

outcome model, even if we completely misspecified the model

for the time-to-event outcome. In the future, it is interest-

ing to see if we can altogether avoid specifying a model for

the time-to-event outcome, and still make valid and efficient

inference on the secondary response.

Even though we present our method in the setting where

the subcohort is selected as a simple random sample from

the full cohort, our method also works when the subcohort

sampling probabilities depend onZ(e.g., stratified random

sampling with respect toZ). This is because our method only

involves estimating the conditional distributionG

X|Z

. When

the subcohort sampling depends on the survival outcome, our

method cannot be directly applied and needs modification.

This is another area of our future research.

Table 5

Analysis of Sister Study

Proposed ELSRSIPW

Estimate

̂

SDEstimate

̂

SDEstimate

̂

SD

Linear ModelInt0.0200.0610.3940.1580.3370.135

(outcome: BMI)PGE-M0.3610.0520.1940.0570.1860.049

NSAIDS0.3970.0490.2570.1160.3080.100

DRINKING−0.2660.055−0.4340.140−0.4460.120

AGE−0.1110.046−0.2660.115−0.2130.099

Cox ModelPGE-M0.1850.103−0.1830.4480.1900.250

BMI0.0810.1250.2690.4260

.1540.397

NSAIDS0.2490.175−0.4930.7190.5800.666

AGE0.0390.183−0.2270.797−0.1930.601

WAIST0.0370.243−0.0610.979−0.2310.860

SMOKE0.4690.331−17.177040.4491.031

1022Biometrics, September2018

Our approach can be extended in several other ways.

Although we use a linear model (2) for the association between

the secondary responseY

2

and all covariates (X,Z), our

method can be applied to other discrete or continuous sec-

ondary outcomes, as long as the conditional distribution of

Y

2

given (X,Z) is specified parametrically. Our models can

also be easily generalized to the case in which the linear

model (2) and Cox model (3) have different covariates. In

our article, we use the empirical distribution to estimate the

conditional distributionG

X|W

whenWis discrete, and use

the kernel method whenWis continuous. WhenWhas both

discrete and continuous components, the product kernel esti-

mator from Li and Racine (2008) could be used. In addition,

instead of assuming that the error termin the linear model

(2) follows a normal distribution, we are working on a more

robust semiparametric estimation where the only assumption

we make isE(|X,Z)=0.

8. Supplementary Materials

Web Appendices referenced in Section 3 and 5, along with the

sourceRcodes for implementing the proposed methods, are

available with this article at theBiometricswebsite on Wiley

Online Library.

Acknowledgements

The authors thank the editors and referees for very construc-

tive comments that greatly improved the presentation of this

article. This research is part of Yinghao Pan’s Ph.D. disserta-

tion. It is partially supported by grants R01 ES021900 from

the National Institute of Environmental Health Sciences and

P01 CA142538 from the National Cancer Institute.

References

Borgan, O., Langholz, B., Samuelsen, S. O., Goldstein, L., and

Pogoda, J. (2000). Exposure stratified case-cohort designs.

Lifetime Data Analysis6,39

–58.

Chen, K. and Lo, S.-H. (1999). Case-cohort and case-control anal-

ysis with cox’s model.Biometrika86, 755

–764.

Cox, D. R. (1972). Regression models and life tables (with discus-

sion).Journal of the Royal Statistical Society34, 187

–220.

Horvitz, D. G. and Thompson, D. J. (1952). A generalization of

sampling without replacement from a finite universe.Journal

of the American Statistical Association47, 663

–685.

Jiang, Y., Scott, A. J., and Wild, C. J. (2006). Secondary analysis

of case-control data.Statistics in Medicine25, 1323

–1339.

Kang, S. and Cai, J. (2009). Marginal hazards model for case-

cohort studies with multiple disease outcomes.Biometrika

96, 887

–901.

Kim, S., Cai, J., and Lu, W. (2013). More efficient estimators for

case-cohort studies.Biometrika100, 695

–708.

Kim, S., Sandler, D. P., Carswell, G., De Roo, L. A., Parks, C. G.,

Cawthon, R., et al. (2011). Telomere length in peripheral

blood and breast cancer risk in a prospective case-cohort

analysis: results from the sister study.Cancer Causes &

Control22, 1061

–1066.

Kim, S., Taylor, J. A., Milne, G. L., and Sandler, D. P. (2013).

Association between urinary prostaglandin e2 metabolite

and breast cancer risk: A prospective, case–cohort study

of postmenopausal women.Cancer Prevention Research6,

511

–518.

Kulich, M. and Lin, D. (2004). Improving the efficiency of

relative-risk estimation in case-cohort studies.Journal of

the American Statistical Association99, 832

–844.

Lee, A., McMurchy, L., and Scott, A. (1997). Re-using data from

case-control studies.Statistics in Medicine16, 1377

–1389.

Li, Q. and Racine, J. S. (2008). Nonparametric estimation of

conditional cdf and quantile functions with mixed categor-

ical and continuous data.Journal of Business & Economic

Statistics26, 423

–434.

Lin, D. and Zeng, D. (2009). Proper analysis of secondary phe-

notype data in case-control association studies.Genetic

Epidemiology33, 256

–265.

Ma, Y. and Carroll, R. J. (2016). Semiparametric estimation in the

secondary analysis of case–control studies.Journal of the

Royal Statistical Society, Series B (Statistical Methodology)

78, 127

–151.

Monsees, G. M., Tamimi, R. M., and Kraft, P. (2009). Genome-

wide association scans for secondary traits using case-control

samples.Genetic Epidemiology33, 717

–728.

Morris, P. G., Hudis, C. A., Giri, D., Morrow, M., Falcone, D.

J., Zhou, X. K., et al. (2011). Inflammation and increased

aromatase expression occur in the breast tissue of obese

women with breast cancer.Cancer Prevention Research4,

1021

–1029.

Prentice, R. L. (1986). A case-cohort design for epidemiologic

cohort studies and disease prevention trials.Biometrika73,

1

–11.

Richardson, D. B., Rzehak, P., Klenk, J., and Weiland, S.

K. (2007). Analyses of case–control data for additional

outcomes.Epidemiology18, 441

–445.

Robins, J. M., Rotnitzky, A., and Zhao, L. P. (1994). Estimation of

regression coefficients when some regressors are not always

observed.Journal of the American Statistical Association

89, 846

–866.

Saarela, O., Kulathinal, S., Arjas, E., and L ̈a ̈ar ̈a, E. (2008).

Nested case–control data utilized for multiple outcomes: A

likelihood approach and alternatives.Statistics in Medicine

27, 5991

–6008.

Saarela, O., Kulathinal, S., and Karvanen, J. (2012). Secondary

analysis under cohort sampling designs using conditional

likelihood.Journal of Probability and Statistics2012,

Article ID 931416.

Salim, A., Xiangmei, M., Jialiang, L., and Reilly, M. (2014). A

maximum likelihood method for secondary analysis of nested

case-control data.Statistics in Medicine33, 1842

–1852.

Self, S. G. and Prentice, R. L. (1988). Asymptotic distribution

theory and efficiency results for case-cohort studies.The

Annals of Statistics16,64

–81.

Subbaramaiah, K., Howe, L. R., Bhardwaj, P., Du, B., Gravaghi,

C., Yantiss, R. K., et al. (2011). Obesity is associated with

inflammation and elevated aromatase expression in the

mouse mammary gland.Cancer Prevention Research4,

329

–346.

Wang, X. and Zhou, H. (2006). A semiparametric empirical like-

lihood method for biased sampling schemes with auxiliary

covariates.Biometrics62, 1149

–1160.

Wang, X. and Zhou, H. (2010). Design and inference for cancer

biomarker study with an outcome and auxiliary-dependent

subsampling.Biometrics66, 502

–511.

Wei, J., Carroll, R. J., M ̈uller, U. U., Keilegom, I. V., and

Chatterjee, N. (2013). Robust estimation for homoscedastic

regression in the secondary analysis of case–control data.

Journal of the Royal Statistical Society, Series B (Statistical

Methodology)75, 185

–206.

Received January2017. Revised November2017.

Accepted November2017.