# Secondary Analysis under Cohort Sampling Designs Using Conditional Likelihood

Authors: Olli Saarela; Sangita Kulathinal; Juha Karvanen
Journal: Journal of Probability and Statistics
 (2012)

---

## Abstract

Under cohort sampling designs, additional covariate data are collected on cases of a specific type and a randomly selected subset of noncases, primarily for the purpose of studying associations with a time-to-event response of interest. With such data available, an interest may arise to reuse them for studying associations between the additional covariate data and a secondary non-time-to-event response variable, usually collected for the whole study cohort at the outset of the study. Following earlier literature, we refer to such a situation as secondary analysis. We outline a general conditional likelihood approach for secondary analysis under cohort sampling designs and discuss the specific situations of case-cohort and nested case-control designs. We also review alternative methods based on full likelihood and inverse probability weighting. We compare the alternative methods for secondary analysis in two simulated settings and apply them in a real-data example.


---

## Body

## 1. Introduction

Cohort sampling designs are two-phase epidemiological study designs where information on time-to-event outcomes of interest over a followup period and some basic covariate data are collected on the whole first-phase study group, referred to as a cohort, and in the second phase, more expensive or difficult-to-obtain additional covariate data are collected only on a subset of the study cohort. This usually comprises the cases, that is, individuals with a disease event of interest during the followup, and a randomly selected subset of noncases. Examples are the case-cohort [1–3] and nested case-control [4, 5] designs. Primarily, such designs are applied for the purpose of studying associations between the time-to-event outcomes and the covariates collected in the second phase. However, with such data having been collected, an interest frequently arises to reuse it for studying associations between the second-phase covariates and the other available covariate data. For instance, the covariates collected in the second phase could be genotypes, while the other covariates may be various phenotype measurements carried out at the outset of the followup period for the whole cohort. The interest would then be to explain a phenotypic response with the genetic covariates. Following Jiang et al. [6] and Lin and Zeng [7], we refer to such a situation as secondary analysis. Here, we concentrate specifically on non-time-to-event secondary outcomes. Analysis of secondary time-to-event outcomes under the nested case-control design has been considered previously by Saarela et al. [8] and Salim et al. [9].As our motivating example, we consider here a single cohort which was used in a larger meta-analysis of association between the European lactase persistence genotype and body mass index (BMI) [10], the latter being a secondary outcome in the cohort study in question. The cohort consists of 5073 men aged 55–77 years from southern and western Finland, who originally formed the placebo group of the ATBC cancer prevention study [11]. Whole blood samples of the participants were taken between 1992 and 1993, which is here considered as the baseline of the cohort, with followup for cardiovascular disease events and all-cause mortality available until the end of year 1999. There is no loss to followup, so the only censoring present is of type I due to end of the followup period. This cohort is a part of MORGAM project, an international pooling of cardiovascular cohorts [12]. Genotype data (including the lactase persistence SNP rs4988235) under this project have been collected under a case-cohort design described in detail by [13] and herein in Section 4.3.1. Given such data, our aim is to estimate the association between the lactase persistence genotype and BMI making use of genotype data collected on both the random subcohort and cases of all-cause mortality.Secondary analysis of case-control data has been studied previously, using profile likelihood [14], inverse selection probability weighting methods [15–17], or retrospective likelihood [6, 7]. However, to the best of our knowledge, a systematic discussion on secondary analysis under cohort sampling designs has been lacking, which we will aim to rectify here by discussing alternative approaches for such an analysis under a generic two-phase study design. We will briefly review the full likelihood approach which utilizes all observed data (Section 2), as well as pseudolikelihoods based on inverse selection probability weighting (Section 3). For these approaches, we propose a conditional likelihood-based alternative (Section 4), restricted to the fully observed second-phase study group. Conditional likelihood inference under cohort sampling designs has been studied previously for the analysis of the primary time-to-event outcome by Langholz and Goldstein [18] and Saarela and Kulathinal [19]; here, we extend these methods to the secondary analysis setting. The main interest is in continuous secondary outcomes, though the approach would also be valid for categorical responses. As special cases of the general setting, we consider case-cohort and nested case-control designs. As extensions to the basic setting, we consider treatment of missing second-phase covariate data and adjustment for left truncation in the case of incident time-to-event outcomes (Section 5). In Section 6, we present two simulation studies, first comparing the efficiencies of the alternative approaches and then demonstrating the potential adverse effects of small sampling fraction in full likelihood inference. We also carry out the analysis in the real-data example using all three alternative methods. As the model for the continuous secondary response variable, in addition to the customary normal distribution, we consider more flexible model specifications, thus aiming to incorporate residual-based model fit diagnostics into the model itself. While other generalizations for the normal distribution have been proposed, we adopt here the four-parameter normal-polynomial quantile mixture [20], which includes the normal distribution as a special case.


## 2. Notation, Assumptions, and Full Likelihood

To cover our motivating example and also the general case, we introduce first some notation. Let the set𝒞≡{1,…,N} represent the individuals in the cohort. Primary outcome in the study is a time-to-event outcome characterized by random variables (Ti,Ei), i∈𝒞, where Ti corresponds to event time and Ei to the event type of individual i, with Ei=0 indicating a censoring event and Ei=1 a death due to any cause. Extension to incident nonfatal outcomes and multiple outcome types is considered separately in Section 5.2. A secondary outcome of interest Yi, here BMI, is observed on all study participants at the outset of the study. In addition, there may be other covariates Xi, available on all i∈𝒞, relevant to be included in the analysis. In the present example, these comprise only the age at the start of the followup. Additional covariate data (here the lactase persistence genotype) Zi are collected only on the second-phase study group 𝒪≡{i:Ri=1}⊆𝒞, specified by the inclusion indicators Ri∈{0,1}, analogously to the survey response/nonresponse setting of Rubin [21]. We will henceforth use vector notations of the type Z𝒪≡{Zi:i∈𝒪} to represent data obtained on different subsets of individuals. We will not make a distinction between random variables and their realized values in the notation when this is clear from the context. The observed data as a whole are then represented as (R𝒞,T𝒞,E𝒞,X𝒞,Y𝒞,Z𝒪). We are interested in the model P(Yi∣Xi,Zi,β) for the secondary outcome, more precisely, the parameter describing the association between Yi and Zi, which in our case correspond to BMI and the genotypic covariate.We assume that the first-phase sampling mechanism has been unconfounded in the sense of Rubin [21, page 36], so that we may ignore the first-phase sampling mechanism in all subsequent analyses. This means that the cohort recruitment (possibly through survey sampling) and possible nonresponse depend only on the Xi-covariates; then if all subsequent analyses are conditional on Xi, the selection mechanism may be ignored. In contrast, the second phase sampling may be outcome dependent; the second-phase sampling mechanism is specified by the joint probability distribution for the set of indicator variables R𝒞. This is assumed to be unconfounded with Z𝒞, that is, R𝒞⊥(Z𝒞,θ)∣T𝒞,E𝒞,X𝒞,Y𝒞. What follows could be further generalized by assuming the sampling mechanism to be ignorable so that R𝒞⊥(Z𝒞∖𝒪,θ)∣T𝒞,E𝒞,X𝒞,Y𝒞,Z𝒪, but since most common cohort sampling mechanisms go under the former assumption, and it will simplify the exposition, we will proceed with that. In addition, we assume the random vectors (Ti,Ei,Yi,Xi,Zi) either to be (infinitely) exchangeable over unit indices i (cf. [21, page 40]), or equivalently, conditionally independent given the collection of all relevant parameters θ. It should be noted that the exchangeability assumption needs not to be extended to the inclusion indicators Ri. Now, following, for example, Saarela et al. [8], we can write a full (observed data) likelihood expression(2.1)P(RC,TC,EC,YC,ZO∣XC,θ)∝(α,β,γ)∏i∈OP(Ti,Ei∣Xi,Yi,Zi,α)P(Yi∣Xi,Zi,β)P(Zi∣Xi,γ)×∏i∈C∖O∫ziP(Ti,Ei∣Xi,Yi,zi,α)P(Yi∣Xi,zi,β)P(dzi∣Xi,γ).
In the secondary analysis situation, only the parameters β are of interest, while α and γ remain present as nuisance parameters. There are potential drawbacks related to the above likelihood expression; it requires integration over the unobserved covariate data in the set 𝒞∖𝒪 (which may be large compared to 𝒪), as well as modeling of the population distribution of Zi. If possible, we would like to avoid this due to the required estimation of the nuisance parameters γ and the risk of model misspecification. Furthermore, observed data likelihoods may become sensitive to misspecification of the model for the response variable; the missing data can act similarly to extra parameters, and the actual model parameters may lose their intended interpretation. This is a real problem especially in cohort sampling designs with a rare event of interest, since the proportion of uncollected covariate data in the study cohort may then be very high. We demonstrate such a situation with simplified simulation example in Section 6.2.


## 3. Methods Based on Inverse Probability Weighting

Valid estimates for the parameters of the secondary outcome model could alternatively be obtained by using inverses of the first-order inclusion probabilitiesP(Ri=1∣T𝒞,E𝒞,X𝒞,Y𝒞) (assumed here to be known) as weights (e.g., [22, 23]). The weighted log-likelihood function, approximating the corresponding complete data log-likelihood, to be used for the secondary analysis is(3.1)∑i∈Olog⁡P(Yi∣Xi,Zi,β)P(Ri=1∣TC,EC,XC,YC)=∑i∈C1{Ri=1}log⁡P(Yi∣Xi,Zi,β)P(Ri=1∣TC,EC,XC,YC).
In a typical cohort sampling design, all cases would receive unit weights, while noncases selected to the set 𝒪 would receive weights greater than one, inverse to their probability to be included in the set of controls/subcohort. While this kind of weighting results in unbiased estimation of the parameters of interest, it will potentially result in reduced efficiency compared to the full and conditional likelihood approaches, since in (3.1) cases receive smaller weights compared to noncases irrespective of whether the case status is actually associated to the secondary outcome or the covariate of interest, and thus the information in these observations may not be fully utilized. We will demonstrate this in a simulation study in Section 6.1. Theoretical justification for estimating function (3.1), as well as variance estimation is discussed in Appendix A. A variation of the above Horvitz-Thompson [24] type of weighting would be poststratification-based estimation (e.g., [15, 25]), with the stratification carried out over the relevant first-phase variables.


## 4. Conditional Likelihood Inference under Cohort Sampling Designs

### 4.1. Definition

A very general definition for conditional likelihood is given by Cox and Hinkley [26, pages 16-17] and Cox [27, page 269] as follows. Let U be a random vector corresponding to all observed data, and let this be partitioned into relevant subsets  U=(V,W), the transformation not depending on unknown parameters θ. The conditional likelihood given the realization V=v is then the conditional probability distribution P(W∈dw∣V=v,θ) with respect to θ. Few generally accepted rules for choosing the partitioning can be found from the literature, one of these being conditioning on an ancillary statistic or something close to that in order to eliminate nuisance parameters (e.g., [28, 29]). In contrast, although working under the same general definition, here we condition on a sampling mechanism, in order to restrict the analysis into a subgroup for which a useful likelihood expression can be written. Generally such a conditioning may lose information on θ, but will nevertheless produce valid estimates. The rules we set out for choosing the partitioning to construct a conditional likelihood under the two-phase study setting are as follows.   (1)
Condition on the sampling mechanism, that is, the set of inclusion indicatorsR𝒞, which produced the second-phase study group 𝒪 onto which the analysis is to be restricted. Do not condition on any further information on the sampling mechanism, such as inclusion in the subcohort or the set of controls, since this can easily lead to overconditioning which will lose a lot of information. In our notation, such additional information on the sampling mechanism is implicitly included in W and, given the assumptions stated in Section 2, will cancel out of the resulting likelihood expression.(2)
Other observed variables may be placed intoV or W at will, depending on the parameters of interest. For instance, if the parameters γ are not of interest, we may condition on Z𝒪 by placing it in V. If the parameters γ need to be estimated, Z𝒪 may be placed in W.(3)
We must haveP(W∈dw∣V=v,θ)P(V∈dv∣θ)=P(U∈du∣θ), that is, all the relevant observed variables must either be modeled or conditioned upon.Applying these conditioning rules will reproduce the conditional likelihood expressions obtained previously in special cases of the current framework by Langholz and Goldstein [18] and Saarela and Kulathinal [19]. The proposed conditional likelihood framework also includes the familiar retrospective likelihood, often suggested for analysis under case-control designs (e.g., [30, 31]), as a special case (see Appendix B). Whereas Langholz and Goldstein [18] considered only the special case of logistic likelihoods, here we aim to first derive the likelihood expressions in the general case before substituting in any specific models.


### 4.2. Conditional Likelihood Expression: The General Case

Following the stated rules, and making the same general assumptions as in Section2, we proceed by partitioning the observed data into relevant subsets W≡(T𝒪,E𝒪,Y𝒪) and V≡(R𝒞,T𝒞∖𝒪,E𝒞∖𝒪,Y𝒞∖𝒪,X𝒞,Z𝒪) and, using a shorthand notation of Qi≡(Ti,Ei,Yi) for all outcome variables, work with a conditional likelihood   (4.1)P(QO∣RC,QC∖O,XC,ZO,θ)=P(RC∣QC,XC,ZO,θ)P(QC,ZO∣XC,θ)P(RC∣QC∖O,XC,ZO,θ)P(QC∖O,ZO∣XC,θ)∝θP(QC,ZO∣XC,θ)P(RC∣QC∖O,XC,ZO,θ)P(QC∖O,ZO∣XC,θ),
where the proportionality follows from the assumption on unconfounded second-phase sampling mechanism. In Appendix C, we show that such a conditional likelihood indeed has the properties of a likelihood, that is, a score function with zero mean and variance equal to the Fisher information. Asymptotic normality of the maximum likelihood estimators based on the conditional likelihood follows obviously from these results in the special case of Bernoulli sampling (Section 4.3.1), although in the general case, this may require further assumptions on the sampling mechanism, such as asymptotic independence of the inclusion indicators Ri.The ratio of the numerator and the second term in the denominator can be further written as(4.2)P(QC,ZO∣XC,θ)P(QC∖O,ZO∣XC,θ)=∏i∈OP(Qi,Zi∣Xi,θ)∏i∈C∖O∫ziP(Qi,Zi∈dzi∣Xi,θ)∏i∈O∫qiP(Qi∈dqi,Zi∣Xi,θ)∏i∈C∖O∫ziP(Qi,Zi∈dzi∣Xi,θ)=∏i∈OP(Qi∣Xi,Zi,θ).
Here, the product forms follow from the exchangeability assumption for the random vectors (Qi,Xi,Zi). Thus, we have obtained   (4.3)P(QO∣RC,QC∖O,XC,ZO,θ)=∏i∈OP(Qi∣Xi,Zi,θ)P(RC∣QC∖O,XC,ZO,θ),
where the numerator factors into the parametric models P(Qi∣Xi,Zi,θ)=P(Ti,Ei∣Xi,Yi,Zi,α)P(Yi∣Xi,Zi,β). The denominator is the conditional likelihood correction term (sometimes called ascertainment correction, e.g., Ma et al. [32]). Its specific form depends on the second-phase sampling mechanism, and in the general case, it does not reduce into a product form. It depends on the model parameters, since it is not conditioned on all of (T𝒞,E𝒞,Y𝒞,X𝒞). Generally, the challenge in conditional likelihood correction terms is in representing them in terms of parameters estimable from the data. Here, the term can be written as   (4.4)P(RC∣QC∖O,XC,ZO,θ)=∫qOP(RC∣QC,XC)P(QO∈dqO∣QC∖O,XC,ZO,θ),
where P(Q𝒪∣Q𝒞∖𝒪,X𝒞,Z𝒪,θ) is given by (4.2), and we have   (4.5)P(RC∣QC∖O,XC,ZO,θ)=∫yi:i∈O∫ti:i∈O∑ei:i∈OP(RC∣TC,EC,YC,XC)×∏i∈O[P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)].
Here, the first term inside the integral specifies the sampling mechanism and is assumed to be known. The obtained likelihood expression utilizes all observed data on covariates Zi, and it can be easily seen that if there is no association between (Ti,Ei) and Zi and between Zi and Xi, we lose no information relevant to learning about the association between Yi and Zi, when the conditional likelihood (4.3) is used instead of the corresponding full likelihood (2.1). This was also demonstrated in the simulation study of Saarela and Kulathinal [19]. On the other hand, if the other observed data do give significant information on the unobserved covariate values Z𝒞∖𝒪, efficiency could potentially be improved by using the full likelihood which utilizes all observed data, with the cost of having to specify a model for Zi. In Section 6.1, the efficiencies of expressions (2.1) and (4.3) are compared in a simulated setting.


### 4.3. Special Cases: Case-Cohort and Nested Case-Control Designs

#### 4.3.1. Case-Cohort/Bernoulli Sampling

Here, we are mainly interested in a variation of the “efficient case-cohort design” suggested by Kim and De Gruttola [33, pages 155-156] as an alternative to sampling within strata. To improve efficiency in sampling of the subcohort, here the distribution of some key covariates (in the present notation either X𝒞 or both (X𝒞,Y𝒞)) in the subcohort is approximately matched to that of the cases by first fitting a logistic regression model, say, logit{P(Ei=1∣Xi,μ)}=μ′Xi, and then selecting the subcohort using Bernoulli sampling with the probabilities πi≡1/(1+exp{-μ̂′Xi}), independently of the case status. We then have P(Ri∣T𝒞,E𝒞,Y𝒞,X𝒞,μ̂)=P(Ri∣Ti,Ei,Yi,Xi,μ̂), where in practice, we make the approximation P(Ri∣Ti,Ei,Yi,Xi,μ̂)≈P(Ri∣Ti,Ei,Yi,Xi,μ) and will subsequently suppress the sampling mechanism parameters from the notation. The selection probabilities may also be rescaled to give a desired expected subcohort size m as πi*≡mπi/∑i∈𝒞πi≈mπi/(NP(Ei=1)). More generally, the subcohort selection probability πi may be any known function of (Ti,Ei,Yi,Xi), that is, πi≡π(Ti,Ei,Yi,Xi), with the sampling design specified by P(R𝒞∣T𝒞,E𝒞,Y𝒞,X𝒞)=∏i∈𝒞P(Ri∣Ti,Ei,Yi,Xi)=∏i∈𝒞π(Ti,Ei,Yi,Xi)Ri(1-π(Ti,Ei,Yi,Xi))1-Ri. The special case of Bernoulli sampling is of interest, because here the product of the first-order inclusion probabilities specifies the joint distribution of the inclusion indicators (see, e.g., [34, pages 62-63]), which considerably simplifies the analysis using conditional likelihood. In the case of stratified without replacement sampling, the following could only be interpreted as a first-order approximation when the sampling fractions are small. In the Bernoulli sampling case, the conditional likelihood correction term reduces into the product form   (4.6)P(RC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)=∏i∈C∖OP(Ri=0∣Ti,Ei,Yi,Xi)×∫yi:i∈O∫ti:i∈O∑ei:i∈O∏i∈O[P(Ri=1∣ti,ei,yi,Xi)×P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)]∝(α,β)∏i∈O∫yi∫ti∑ei[P(Ri=1∣ti,ei,yi,Xi)×P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)].
Let ℰ≡{i:Ei=1} denote the set of individuals who died during the followup, and let d≡|ℰ|. In the case-cohort design discussed in Kulathinal and Arjas [35], Kulathinal et al. [13], and Saarela and Kulathinal [19], all cases are selected to the case-cohort set, so that P(Ri=1∣Ti,Ei=1,Yi,Xi)=1, i∈ℰ, while to increase the efficiency of the design, the subcohort is selected with probabilities which depend on the age bi at the start of the followup (included in the Xi covariates) through a logistic model as discussed above, giving the inclusion probability for non-cases as P(Ri=1∣Ti,Ei=0,Yi,Xi)=π(bi), i∈𝒞∖ℰ. Under Bernoulli sampling, the realized sample size n≡|𝒪| is random, with the expected sample size in the present example given by E(n∣E𝒞,X𝒞)=d+∑i∈𝒞∖ℰπ(bi). In the case of a mortality outcome and type I censoring at predetermined times ci, with the observed time given by Ti≡min(T̃i,ci), where T̃i is the latent time of death, (4.6) further simplifies into (see also [19, pages 12-13])(4.7)P(RC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)∝(α,β)∏i∈O∫yi∫ti∈[0,ci][P(Ri=1∣ti,Ei=1,yi,Xi)P(dti,Ei=1∣Xi,yi,Zi,α)+P(Ri=1∣ti,Ei=0,yi,Xi)P(dti,Ei=0∣Xi,yi,Zi,α)]×P(dyi∣Xi,Zi,β)=∏i∈O∫yi[P(0≤Ti<ci,Ei=1∣Xi,yi,Zi,α)+π(bi)P(Ti∈dci,Ei=0∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β)=∏i∈O∫yi[P(0≤T̃i<ci∣Xi,yi,Zi,α)+π(bi)P(T̃i≥ci∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β)=∏i∈O∫yi[1-(1-π(bi))P(T̃i≥ci∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β).
Hence, (4.6) could be represented in terms of the survival function. Suppose now that in addition to the type I censoring due to the end of the followup period, there may be random censoring during the followup; for instance, this may be the case if the outcome in the case-cohort design (Ei=1) is not all-cause mortality, but, say, mortality due to cardiovascular diseases. Deaths due to other causes, denoted as Ei=2, then appear as censoring. As before, let T̃i denote the latent time of death, with Ei indicating the type of death or end of followup period. Similarly as above, we get   (4.8)P(RC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)=∏i∈O∫yi[P(0≤Ti<ci,Ei=1∣Xi,yi,Zi,α)+π(bi)P(0≤Ti<ci,Ei=2∣Xi,yi,Zi,α)+π(bi)P(Ti∈dci,Ei=0∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β)=∏i∈O∫yi[P(0≤T̃i<ci,Ei=1∣Xi,yi,Zi,α)+π(bi)P(0≤T̃i<ci,Ei=2∣Xi,yi,Zi,α)+π(bi)P(T̃i≥ci∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β).
Here, the probabilities P(0≤T̃i<ci,Ei=k∣Xi,Yi,Zi,α), k∈{1,2}, are given by the cumulative incidence functions and P(T̃i≥ci∣Xi,Yi,Zi,α) by the joint survival function of a competing risks survival model for the two types of death (e.g., [36, pages 251-252]). These are in principle identifiable from the observed data for the set 𝒪, since the subcohort has been selected independently of the case status and thus will include a number of other deaths. On the other hand, if this number is very small, the middle term in the above sum contributes little to the correction term, and thus unstable estimation of the corresponding parameters does not necessarily hinder the estimation of the parameters of interest.


#### 4.3.2. Risk Set Sampling

Consider now a nested case-control sampling mechanism in which all cases are selected to the case-control set with probability one, andm controls per case j are selected using without replacement sampling from the risk set ℛj≡{i:Ai(Tj)=1}∖{j} (of size rj≡|ℛj|=∑i∈𝒞Ai(Tj)-1), where Ai(Tj) is the at-risk indicator for subject i at event time Tj. This is carried out independently for all cases j∈ℰ. Let the sampled set of time-matched controls for case j be denoted as 𝒮j⊆ℛj, some of which may also be future cases since the sampling is carried out without regard to the future case status of the individuals in the risk set. Let 𝒮ℰ≡{𝒮j:j∈ℰ} denote the collection of the sampled risk sets and 𝒮≡⋃j∈ℰ𝒮j the pooled set of sampled controls. Noting that knowing all of (T𝒞,E𝒞) specifies the risk sets ℛℰ≡{ℛj:j∈ℰ} as well as the order of the events, the risk set sampling mechanism is specified by the joint probability distribution   (4.9)P(RC∣TC,EC,YC,XC)=P(RC∣RE,E)=∑SE:|Sj|=m,j∈EP(RC∣RE,SE,E)P(SE∣RE,E)=∑SE:|Sj|=m,j∈E1{S∪E=O}∏j∈EP(Sj∣Rj)=∏j∈E1(rjm)∑SE:|Sj|=m,j∈E1{S∪E=O}.
If the sampling has been carried out within strata specified by the covariates Y𝒞 and/or X𝒞 through some function g(Xi,Yi), the risk sets would be defined within strata as ℛj={i:Ai(Tj)=1,g(Xi,Yi)=g(Xj,Yj)}∖{j}, and the above reasoning would apply similarly with the redefined risk sets. The conditional likelihood correction term now becomes   (4.10)P(RC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)=∫yi:i∈O∫ti:i∈O∑ei:i∈O∏j∈E1(rjm)∑SE:|Sj|=m,j∈E1{S∪E=O}×∏i∈O[P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)],
which does not reduce into a product form. Exact numerical evaluation of this term would require enumeration of all possible combinations of sampled risk sets of size m which would result in the observed case-control set 𝒪. This is unlikely to be feasible in practice with realistic sample sizes, which is why consideration of approximations may be necessary. Samuelsen [23] showed the inclusion indicators under risk set sampling to be asymptotically pairwise uncorrelated, with the first-order inclusion probabilities for the noncases i∈𝒞∖ℰ given by P(Ri=1∣T𝒞,E𝒞∖{i},Ei=0)=1-∏j∈ℰ[1-mAi(Tj)/rj]. Thus, it might be tempting to replace (4.9) with ∏i∈𝒞P(Ri∣T𝒞,E𝒞,Y𝒞,X𝒞). However, the properties of such an approximation will require further research.


## 4.1. Definition

A very general definition for conditional likelihood is given by Cox and Hinkley [26, pages 16-17] and Cox [27, page 269] as follows. Let U be a random vector corresponding to all observed data, and let this be partitioned into relevant subsets  U=(V,W), the transformation not depending on unknown parameters θ. The conditional likelihood given the realization V=v is then the conditional probability distribution P(W∈dw∣V=v,θ) with respect to θ. Few generally accepted rules for choosing the partitioning can be found from the literature, one of these being conditioning on an ancillary statistic or something close to that in order to eliminate nuisance parameters (e.g., [28, 29]). In contrast, although working under the same general definition, here we condition on a sampling mechanism, in order to restrict the analysis into a subgroup for which a useful likelihood expression can be written. Generally such a conditioning may lose information on θ, but will nevertheless produce valid estimates. The rules we set out for choosing the partitioning to construct a conditional likelihood under the two-phase study setting are as follows.   (1)
Condition on the sampling mechanism, that is, the set of inclusion indicatorsR𝒞, which produced the second-phase study group 𝒪 onto which the analysis is to be restricted. Do not condition on any further information on the sampling mechanism, such as inclusion in the subcohort or the set of controls, since this can easily lead to overconditioning which will lose a lot of information. In our notation, such additional information on the sampling mechanism is implicitly included in W and, given the assumptions stated in Section 2, will cancel out of the resulting likelihood expression.(2)
Other observed variables may be placed intoV or W at will, depending on the parameters of interest. For instance, if the parameters γ are not of interest, we may condition on Z𝒪 by placing it in V. If the parameters γ need to be estimated, Z𝒪 may be placed in W.(3)
We must haveP(W∈dw∣V=v,θ)P(V∈dv∣θ)=P(U∈du∣θ), that is, all the relevant observed variables must either be modeled or conditioned upon.Applying these conditioning rules will reproduce the conditional likelihood expressions obtained previously in special cases of the current framework by Langholz and Goldstein [18] and Saarela and Kulathinal [19]. The proposed conditional likelihood framework also includes the familiar retrospective likelihood, often suggested for analysis under case-control designs (e.g., [30, 31]), as a special case (see Appendix B). Whereas Langholz and Goldstein [18] considered only the special case of logistic likelihoods, here we aim to first derive the likelihood expressions in the general case before substituting in any specific models.


## 4.2. Conditional Likelihood Expression: The General Case

Following the stated rules, and making the same general assumptions as in Section2, we proceed by partitioning the observed data into relevant subsets W≡(T𝒪,E𝒪,Y𝒪) and V≡(R𝒞,T𝒞∖𝒪,E𝒞∖𝒪,Y𝒞∖𝒪,X𝒞,Z𝒪) and, using a shorthand notation of Qi≡(Ti,Ei,Yi) for all outcome variables, work with a conditional likelihood   (4.1)P(QO∣RC,QC∖O,XC,ZO,θ)=P(RC∣QC,XC,ZO,θ)P(QC,ZO∣XC,θ)P(RC∣QC∖O,XC,ZO,θ)P(QC∖O,ZO∣XC,θ)∝θP(QC,ZO∣XC,θ)P(RC∣QC∖O,XC,ZO,θ)P(QC∖O,ZO∣XC,θ),
where the proportionality follows from the assumption on unconfounded second-phase sampling mechanism. In Appendix C, we show that such a conditional likelihood indeed has the properties of a likelihood, that is, a score function with zero mean and variance equal to the Fisher information. Asymptotic normality of the maximum likelihood estimators based on the conditional likelihood follows obviously from these results in the special case of Bernoulli sampling (Section 4.3.1), although in the general case, this may require further assumptions on the sampling mechanism, such as asymptotic independence of the inclusion indicators Ri.The ratio of the numerator and the second term in the denominator can be further written as(4.2)P(QC,ZO∣XC,θ)P(QC∖O,ZO∣XC,θ)=∏i∈OP(Qi,Zi∣Xi,θ)∏i∈C∖O∫ziP(Qi,Zi∈dzi∣Xi,θ)∏i∈O∫qiP(Qi∈dqi,Zi∣Xi,θ)∏i∈C∖O∫ziP(Qi,Zi∈dzi∣Xi,θ)=∏i∈OP(Qi∣Xi,Zi,θ).
Here, the product forms follow from the exchangeability assumption for the random vectors (Qi,Xi,Zi). Thus, we have obtained   (4.3)P(QO∣RC,QC∖O,XC,ZO,θ)=∏i∈OP(Qi∣Xi,Zi,θ)P(RC∣QC∖O,XC,ZO,θ),
where the numerator factors into the parametric models P(Qi∣Xi,Zi,θ)=P(Ti,Ei∣Xi,Yi,Zi,α)P(Yi∣Xi,Zi,β). The denominator is the conditional likelihood correction term (sometimes called ascertainment correction, e.g., Ma et al. [32]). Its specific form depends on the second-phase sampling mechanism, and in the general case, it does not reduce into a product form. It depends on the model parameters, since it is not conditioned on all of (T𝒞,E𝒞,Y𝒞,X𝒞). Generally, the challenge in conditional likelihood correction terms is in representing them in terms of parameters estimable from the data. Here, the term can be written as   (4.4)P(RC∣QC∖O,XC,ZO,θ)=∫qOP(RC∣QC,XC)P(QO∈dqO∣QC∖O,XC,ZO,θ),
where P(Q𝒪∣Q𝒞∖𝒪,X𝒞,Z𝒪,θ) is given by (4.2), and we have   (4.5)P(RC∣QC∖O,XC,ZO,θ)=∫yi:i∈O∫ti:i∈O∑ei:i∈OP(RC∣TC,EC,YC,XC)×∏i∈O[P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)].
Here, the first term inside the integral specifies the sampling mechanism and is assumed to be known. The obtained likelihood expression utilizes all observed data on covariates Zi, and it can be easily seen that if there is no association between (Ti,Ei) and Zi and between Zi and Xi, we lose no information relevant to learning about the association between Yi and Zi, when the conditional likelihood (4.3) is used instead of the corresponding full likelihood (2.1). This was also demonstrated in the simulation study of Saarela and Kulathinal [19]. On the other hand, if the other observed data do give significant information on the unobserved covariate values Z𝒞∖𝒪, efficiency could potentially be improved by using the full likelihood which utilizes all observed data, with the cost of having to specify a model for Zi. In Section 6.1, the efficiencies of expressions (2.1) and (4.3) are compared in a simulated setting.


## 4.3. Special Cases: Case-Cohort and Nested Case-Control Designs

### 4.3.1. Case-Cohort/Bernoulli Sampling

Here, we are mainly interested in a variation of the “efficient case-cohort design” suggested by Kim and De Gruttola [33, pages 155-156] as an alternative to sampling within strata. To improve efficiency in sampling of the subcohort, here the distribution of some key covariates (in the present notation either X𝒞 or both (X𝒞,Y𝒞)) in the subcohort is approximately matched to that of the cases by first fitting a logistic regression model, say, logit{P(Ei=1∣Xi,μ)}=μ′Xi, and then selecting the subcohort using Bernoulli sampling with the probabilities πi≡1/(1+exp{-μ̂′Xi}), independently of the case status. We then have P(Ri∣T𝒞,E𝒞,Y𝒞,X𝒞,μ̂)=P(Ri∣Ti,Ei,Yi,Xi,μ̂), where in practice, we make the approximation P(Ri∣Ti,Ei,Yi,Xi,μ̂)≈P(Ri∣Ti,Ei,Yi,Xi,μ) and will subsequently suppress the sampling mechanism parameters from the notation. The selection probabilities may also be rescaled to give a desired expected subcohort size m as πi*≡mπi/∑i∈𝒞πi≈mπi/(NP(Ei=1)). More generally, the subcohort selection probability πi may be any known function of (Ti,Ei,Yi,Xi), that is, πi≡π(Ti,Ei,Yi,Xi), with the sampling design specified by P(R𝒞∣T𝒞,E𝒞,Y𝒞,X𝒞)=∏i∈𝒞P(Ri∣Ti,Ei,Yi,Xi)=∏i∈𝒞π(Ti,Ei,Yi,Xi)Ri(1-π(Ti,Ei,Yi,Xi))1-Ri. The special case of Bernoulli sampling is of interest, because here the product of the first-order inclusion probabilities specifies the joint distribution of the inclusion indicators (see, e.g., [34, pages 62-63]), which considerably simplifies the analysis using conditional likelihood. In the case of stratified without replacement sampling, the following could only be interpreted as a first-order approximation when the sampling fractions are small. In the Bernoulli sampling case, the conditional likelihood correction term reduces into the product form   (4.6)P(RC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)=∏i∈C∖OP(Ri=0∣Ti,Ei,Yi,Xi)×∫yi:i∈O∫ti:i∈O∑ei:i∈O∏i∈O[P(Ri=1∣ti,ei,yi,Xi)×P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)]∝(α,β)∏i∈O∫yi∫ti∑ei[P(Ri=1∣ti,ei,yi,Xi)×P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)].
Let ℰ≡{i:Ei=1} denote the set of individuals who died during the followup, and let d≡|ℰ|. In the case-cohort design discussed in Kulathinal and Arjas [35], Kulathinal et al. [13], and Saarela and Kulathinal [19], all cases are selected to the case-cohort set, so that P(Ri=1∣Ti,Ei=1,Yi,Xi)=1, i∈ℰ, while to increase the efficiency of the design, the subcohort is selected with probabilities which depend on the age bi at the start of the followup (included in the Xi covariates) through a logistic model as discussed above, giving the inclusion probability for non-cases as P(Ri=1∣Ti,Ei=0,Yi,Xi)=π(bi), i∈𝒞∖ℰ. Under Bernoulli sampling, the realized sample size n≡|𝒪| is random, with the expected sample size in the present example given by E(n∣E𝒞,X𝒞)=d+∑i∈𝒞∖ℰπ(bi). In the case of a mortality outcome and type I censoring at predetermined times ci, with the observed time given by Ti≡min(T̃i,ci), where T̃i is the latent time of death, (4.6) further simplifies into (see also [19, pages 12-13])(4.7)P(RC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)∝(α,β)∏i∈O∫yi∫ti∈[0,ci][P(Ri=1∣ti,Ei=1,yi,Xi)P(dti,Ei=1∣Xi,yi,Zi,α)+P(Ri=1∣ti,Ei=0,yi,Xi)P(dti,Ei=0∣Xi,yi,Zi,α)]×P(dyi∣Xi,Zi,β)=∏i∈O∫yi[P(0≤Ti<ci,Ei=1∣Xi,yi,Zi,α)+π(bi)P(Ti∈dci,Ei=0∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β)=∏i∈O∫yi[P(0≤T̃i<ci∣Xi,yi,Zi,α)+π(bi)P(T̃i≥ci∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β)=∏i∈O∫yi[1-(1-π(bi))P(T̃i≥ci∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β).
Hence, (4.6) could be represented in terms of the survival function. Suppose now that in addition to the type I censoring due to the end of the followup period, there may be random censoring during the followup; for instance, this may be the case if the outcome in the case-cohort design (Ei=1) is not all-cause mortality, but, say, mortality due to cardiovascular diseases. Deaths due to other causes, denoted as Ei=2, then appear as censoring. As before, let T̃i denote the latent time of death, with Ei indicating the type of death or end of followup period. Similarly as above, we get   (4.8)P(RC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)=∏i∈O∫yi[P(0≤Ti<ci,Ei=1∣Xi,yi,Zi,α)+π(bi)P(0≤Ti<ci,Ei=2∣Xi,yi,Zi,α)+π(bi)P(Ti∈dci,Ei=0∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β)=∏i∈O∫yi[P(0≤T̃i<ci,Ei=1∣Xi,yi,Zi,α)+π(bi)P(0≤T̃i<ci,Ei=2∣Xi,yi,Zi,α)+π(bi)P(T̃i≥ci∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β).
Here, the probabilities P(0≤T̃i<ci,Ei=k∣Xi,Yi,Zi,α), k∈{1,2}, are given by the cumulative incidence functions and P(T̃i≥ci∣Xi,Yi,Zi,α) by the joint survival function of a competing risks survival model for the two types of death (e.g., [36, pages 251-252]). These are in principle identifiable from the observed data for the set 𝒪, since the subcohort has been selected independently of the case status and thus will include a number of other deaths. On the other hand, if this number is very small, the middle term in the above sum contributes little to the correction term, and thus unstable estimation of the corresponding parameters does not necessarily hinder the estimation of the parameters of interest.


### 4.3.2. Risk Set Sampling

Consider now a nested case-control sampling mechanism in which all cases are selected to the case-control set with probability one, andm controls per case j are selected using without replacement sampling from the risk set ℛj≡{i:Ai(Tj)=1}∖{j} (of size rj≡|ℛj|=∑i∈𝒞Ai(Tj)-1), where Ai(Tj) is the at-risk indicator for subject i at event time Tj. This is carried out independently for all cases j∈ℰ. Let the sampled set of time-matched controls for case j be denoted as 𝒮j⊆ℛj, some of which may also be future cases since the sampling is carried out without regard to the future case status of the individuals in the risk set. Let 𝒮ℰ≡{𝒮j:j∈ℰ} denote the collection of the sampled risk sets and 𝒮≡⋃j∈ℰ𝒮j the pooled set of sampled controls. Noting that knowing all of (T𝒞,E𝒞) specifies the risk sets ℛℰ≡{ℛj:j∈ℰ} as well as the order of the events, the risk set sampling mechanism is specified by the joint probability distribution   (4.9)P(RC∣TC,EC,YC,XC)=P(RC∣RE,E)=∑SE:|Sj|=m,j∈EP(RC∣RE,SE,E)P(SE∣RE,E)=∑SE:|Sj|=m,j∈E1{S∪E=O}∏j∈EP(Sj∣Rj)=∏j∈E1(rjm)∑SE:|Sj|=m,j∈E1{S∪E=O}.
If the sampling has been carried out within strata specified by the covariates Y𝒞 and/or X𝒞 through some function g(Xi,Yi), the risk sets would be defined within strata as ℛj={i:Ai(Tj)=1,g(Xi,Yi)=g(Xj,Yj)}∖{j}, and the above reasoning would apply similarly with the redefined risk sets. The conditional likelihood correction term now becomes   (4.10)P(RC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)=∫yi:i∈O∫ti:i∈O∑ei:i∈O∏j∈E1(rjm)∑SE:|Sj|=m,j∈E1{S∪E=O}×∏i∈O[P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)],
which does not reduce into a product form. Exact numerical evaluation of this term would require enumeration of all possible combinations of sampled risk sets of size m which would result in the observed case-control set 𝒪. This is unlikely to be feasible in practice with realistic sample sizes, which is why consideration of approximations may be necessary. Samuelsen [23] showed the inclusion indicators under risk set sampling to be asymptotically pairwise uncorrelated, with the first-order inclusion probabilities for the noncases i∈𝒞∖ℰ given by P(Ri=1∣T𝒞,E𝒞∖{i},Ei=0)=1-∏j∈ℰ[1-mAi(Tj)/rj]. Thus, it might be tempting to replace (4.9) with ∏i∈𝒞P(Ri∣T𝒞,E𝒞,Y𝒞,X𝒞). However, the properties of such an approximation will require further research.


## 4.3.1. Case-Cohort/Bernoulli Sampling

Here, we are mainly interested in a variation of the “efficient case-cohort design” suggested by Kim and De Gruttola [33, pages 155-156] as an alternative to sampling within strata. To improve efficiency in sampling of the subcohort, here the distribution of some key covariates (in the present notation either X𝒞 or both (X𝒞,Y𝒞)) in the subcohort is approximately matched to that of the cases by first fitting a logistic regression model, say, logit{P(Ei=1∣Xi,μ)}=μ′Xi, and then selecting the subcohort using Bernoulli sampling with the probabilities πi≡1/(1+exp{-μ̂′Xi}), independently of the case status. We then have P(Ri∣T𝒞,E𝒞,Y𝒞,X𝒞,μ̂)=P(Ri∣Ti,Ei,Yi,Xi,μ̂), where in practice, we make the approximation P(Ri∣Ti,Ei,Yi,Xi,μ̂)≈P(Ri∣Ti,Ei,Yi,Xi,μ) and will subsequently suppress the sampling mechanism parameters from the notation. The selection probabilities may also be rescaled to give a desired expected subcohort size m as πi*≡mπi/∑i∈𝒞πi≈mπi/(NP(Ei=1)). More generally, the subcohort selection probability πi may be any known function of (Ti,Ei,Yi,Xi), that is, πi≡π(Ti,Ei,Yi,Xi), with the sampling design specified by P(R𝒞∣T𝒞,E𝒞,Y𝒞,X𝒞)=∏i∈𝒞P(Ri∣Ti,Ei,Yi,Xi)=∏i∈𝒞π(Ti,Ei,Yi,Xi)Ri(1-π(Ti,Ei,Yi,Xi))1-Ri. The special case of Bernoulli sampling is of interest, because here the product of the first-order inclusion probabilities specifies the joint distribution of the inclusion indicators (see, e.g., [34, pages 62-63]), which considerably simplifies the analysis using conditional likelihood. In the case of stratified without replacement sampling, the following could only be interpreted as a first-order approximation when the sampling fractions are small. In the Bernoulli sampling case, the conditional likelihood correction term reduces into the product form   (4.6)P(RC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)=∏i∈C∖OP(Ri=0∣Ti,Ei,Yi,Xi)×∫yi:i∈O∫ti:i∈O∑ei:i∈O∏i∈O[P(Ri=1∣ti,ei,yi,Xi)×P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)]∝(α,β)∏i∈O∫yi∫ti∑ei[P(Ri=1∣ti,ei,yi,Xi)×P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)].
Let ℰ≡{i:Ei=1} denote the set of individuals who died during the followup, and let d≡|ℰ|. In the case-cohort design discussed in Kulathinal and Arjas [35], Kulathinal et al. [13], and Saarela and Kulathinal [19], all cases are selected to the case-cohort set, so that P(Ri=1∣Ti,Ei=1,Yi,Xi)=1, i∈ℰ, while to increase the efficiency of the design, the subcohort is selected with probabilities which depend on the age bi at the start of the followup (included in the Xi covariates) through a logistic model as discussed above, giving the inclusion probability for non-cases as P(Ri=1∣Ti,Ei=0,Yi,Xi)=π(bi), i∈𝒞∖ℰ. Under Bernoulli sampling, the realized sample size n≡|𝒪| is random, with the expected sample size in the present example given by E(n∣E𝒞,X𝒞)=d+∑i∈𝒞∖ℰπ(bi). In the case of a mortality outcome and type I censoring at predetermined times ci, with the observed time given by Ti≡min(T̃i,ci), where T̃i is the latent time of death, (4.6) further simplifies into (see also [19, pages 12-13])(4.7)P(RC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)∝(α,β)∏i∈O∫yi∫ti∈[0,ci][P(Ri=1∣ti,Ei=1,yi,Xi)P(dti,Ei=1∣Xi,yi,Zi,α)+P(Ri=1∣ti,Ei=0,yi,Xi)P(dti,Ei=0∣Xi,yi,Zi,α)]×P(dyi∣Xi,Zi,β)=∏i∈O∫yi[P(0≤Ti<ci,Ei=1∣Xi,yi,Zi,α)+π(bi)P(Ti∈dci,Ei=0∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β)=∏i∈O∫yi[P(0≤T̃i<ci∣Xi,yi,Zi,α)+π(bi)P(T̃i≥ci∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β)=∏i∈O∫yi[1-(1-π(bi))P(T̃i≥ci∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β).
Hence, (4.6) could be represented in terms of the survival function. Suppose now that in addition to the type I censoring due to the end of the followup period, there may be random censoring during the followup; for instance, this may be the case if the outcome in the case-cohort design (Ei=1) is not all-cause mortality, but, say, mortality due to cardiovascular diseases. Deaths due to other causes, denoted as Ei=2, then appear as censoring. As before, let T̃i denote the latent time of death, with Ei indicating the type of death or end of followup period. Similarly as above, we get   (4.8)P(RC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)=∏i∈O∫yi[P(0≤Ti<ci,Ei=1∣Xi,yi,Zi,α)+π(bi)P(0≤Ti<ci,Ei=2∣Xi,yi,Zi,α)+π(bi)P(Ti∈dci,Ei=0∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β)=∏i∈O∫yi[P(0≤T̃i<ci,Ei=1∣Xi,yi,Zi,α)+π(bi)P(0≤T̃i<ci,Ei=2∣Xi,yi,Zi,α)+π(bi)P(T̃i≥ci∣Xi,yi,Zi,α)]P(dyi∣Xi,Zi,β).
Here, the probabilities P(0≤T̃i<ci,Ei=k∣Xi,Yi,Zi,α), k∈{1,2}, are given by the cumulative incidence functions and P(T̃i≥ci∣Xi,Yi,Zi,α) by the joint survival function of a competing risks survival model for the two types of death (e.g., [36, pages 251-252]). These are in principle identifiable from the observed data for the set 𝒪, since the subcohort has been selected independently of the case status and thus will include a number of other deaths. On the other hand, if this number is very small, the middle term in the above sum contributes little to the correction term, and thus unstable estimation of the corresponding parameters does not necessarily hinder the estimation of the parameters of interest.


## 4.3.2. Risk Set Sampling

Consider now a nested case-control sampling mechanism in which all cases are selected to the case-control set with probability one, andm controls per case j are selected using without replacement sampling from the risk set ℛj≡{i:Ai(Tj)=1}∖{j} (of size rj≡|ℛj|=∑i∈𝒞Ai(Tj)-1), where Ai(Tj) is the at-risk indicator for subject i at event time Tj. This is carried out independently for all cases j∈ℰ. Let the sampled set of time-matched controls for case j be denoted as 𝒮j⊆ℛj, some of which may also be future cases since the sampling is carried out without regard to the future case status of the individuals in the risk set. Let 𝒮ℰ≡{𝒮j:j∈ℰ} denote the collection of the sampled risk sets and 𝒮≡⋃j∈ℰ𝒮j the pooled set of sampled controls. Noting that knowing all of (T𝒞,E𝒞) specifies the risk sets ℛℰ≡{ℛj:j∈ℰ} as well as the order of the events, the risk set sampling mechanism is specified by the joint probability distribution   (4.9)P(RC∣TC,EC,YC,XC)=P(RC∣RE,E)=∑SE:|Sj|=m,j∈EP(RC∣RE,SE,E)P(SE∣RE,E)=∑SE:|Sj|=m,j∈E1{S∪E=O}∏j∈EP(Sj∣Rj)=∏j∈E1(rjm)∑SE:|Sj|=m,j∈E1{S∪E=O}.
If the sampling has been carried out within strata specified by the covariates Y𝒞 and/or X𝒞 through some function g(Xi,Yi), the risk sets would be defined within strata as ℛj={i:Ai(Tj)=1,g(Xi,Yi)=g(Xj,Yj)}∖{j}, and the above reasoning would apply similarly with the redefined risk sets. The conditional likelihood correction term now becomes   (4.10)P(RC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)=∫yi:i∈O∫ti:i∈O∑ei:i∈O∏j∈E1(rjm)∑SE:|Sj|=m,j∈E1{S∪E=O}×∏i∈O[P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)],
which does not reduce into a product form. Exact numerical evaluation of this term would require enumeration of all possible combinations of sampled risk sets of size m which would result in the observed case-control set 𝒪. This is unlikely to be feasible in practice with realistic sample sizes, which is why consideration of approximations may be necessary. Samuelsen [23] showed the inclusion indicators under risk set sampling to be asymptotically pairwise uncorrelated, with the first-order inclusion probabilities for the noncases i∈𝒞∖ℰ given by P(Ri=1∣T𝒞,E𝒞∖{i},Ei=0)=1-∏j∈ℰ[1-mAi(Tj)/rj]. Thus, it might be tempting to replace (4.9) with ∏i∈𝒞P(Ri∣T𝒞,E𝒞,Y𝒞,X𝒞). However, the properties of such an approximation will require further research.


## 5. Extensions

### 5.1. Missing Second-Phase Covariate Data

Further issue to be considered in practical applications is possible missing covariate data within the set𝒪. For instance, if the covariates to be collected under the cohort sampling design are genotypes, after selection of subjects to be genotyped, it may turn out that the DNA amount or concentration is too low or the genotyping is otherwise unsuccessful. Let M𝒪≡{Mi:i∈𝒪} be a set of indicator variables indicating whether the measurement of Zi turned out to be unsuccessful after selection of subject i into the set 𝒪, and let ℳ≡{i:Mi=1}⊆𝒪. Now Zi is observed only on the set 𝒪∖ℳ. We consider the implications of this separately for all of the three approaches discussed above.


#### 5.1.1. Full Likelihood

If the missingness can be assumed missing at random, that is,M𝒪⊥(Zℳ,θ)∣R𝒞,T𝒞,E𝒞,X𝒞,Y𝒞,Z𝒪∖ℳ, the missingness indicators can again be included in the proportionality constant, and the full likelihood expression becomes (5.1)P(RC,MO,TC,EC,YC,ZO∖M∣XC,θ)∝(α,β,γ)∏i∈O∖MP(Ti,Ei∣Xi,Yi,Zi,α)P(Yi∣Xi,Zi,β)P(Zi∣Xi,γ)×∏i∈C∖O∪M∫ziP(Ti,Ei∣Xi,Yi,zi,α)P(Yi∣Xi,zi,β)P(dzi∣Xi,γ).


#### 5.1.2. Inverse Probability Weighting

If the missingness mechanism would be known (and missing at random), a weighted pseudolikelihood expression for the set𝒪∖ℳ could be obtained by multiplying the original selection probabilities by the probabilities of observing the value of Zi after selection   (5.2)∑i∈O∖Mlog⁡P(Yi∣Xi,Zi,β)P(Mi=0∣RC,TC,EC,XC,YC,ZO∖M)P(Ri=1∣TC,EC,XC,YC).
In practice, however, the missingness mechanism would have to be modeled to obtain estimates for these probabilities.


#### 5.1.3. Conditional Likelihood

By partitioning the observed data intoW=(T𝒪∖ℳ,E𝒪∖ℳ,Y𝒪∖ℳ) and V=(R𝒞,M𝒪,T𝒞∖𝒪∪ℳ,E𝒞∖𝒪∪ℳ,Y𝒞∖𝒪∪ℳ,X𝒞,Z𝒪∖ℳ), we obtain a conditional likelihood (5.3)P(QO∖M∣RC,MO,QC∖O∪M,XC,ZO∖M,θ)∝θP(MO∣RC,QC,XC,ZO∖M,θ)P(MO∣RC,QC∖O∪M,XC,ZO∖M,θ)∏i∈O∖MP(Qi∣Xi,Zi,θ)P(RC∣QC∖O∪M,XC,ZO∖M,θ),
where P(R𝒞∣Q𝒞∖𝒪∪ℳ,X𝒞,Z𝒪∖ℳ,θ) is obtained from (4.5) by replacing the set 𝒪 with 𝒪∖ℳ. Unlike the second-phase sampling mechanism, the missingness mechanism is generally unknown. From (5.3), it can be seen that if we are willing to assume that Zℳ are either missing completely at random or that M𝒪 depends only on X𝒞, the terms involving the missingness indicators cancel out of the likelihood. On the other hand, if the missingness may depend on the response variables Q𝒞, the missingness mechanism would have to be modeled. However, in such a case, it may be easier to work with the partitioning W=(M𝒪,T𝒪,E𝒪,Y𝒪,Z𝒪∖ℳ) and V=(R𝒞,T𝒞∖𝒪,E𝒞∖𝒪,Y𝒞∖𝒪,X𝒞) and model the population distribution of Zi rather than trying to estimate parameters describing the missingness mechanism. This in fact corresponds to what was done by Saarela and Kulathinal [19] and was required there because haplotypes are only partially identifiable from unphased genotype data.


### 5.2. Incident Outcomes and Left Truncation

Previous discussion was specific to a primary mortality outcome using time on study as the main time scale. In this section, we discuss separately how the different methods can accommodate cohort sampling for incident nonfatal primary outcomes. In the analysis of secondary, non-time-to-event outcomes, the presence of left truncation due to exclusion of cases of prevalent disease presents an additional complication. If the parameters of the secondary outcome model correspond to the background population alive at the cohort baseline (rather than to the disease-free population), this additional selection factor requires further adjustment. If the primary outcome is a mortality endpoint, this is not an issue, since then there is no further selection due to prevalent conditions. In likelihood-based adjustment for left truncation, the main time scale of the analysis has to be chosen as age instead of time on study. In survival modeling, it is well known that conditioning on event-free survival until the age at study baseline corresponds to exclusion of the followup time before that (e.g., [37, page 580], and the references therein). However, this is no longer true in the case of missing covariate data (e.g., [8, pages 5997-5998]), or indeed in the analysis of secondary outcomes, as will be demonstrated below. The set notations 𝒞 and 𝒪 are here taken to refer to the disease free cohort and second-phase study group.It should also be noted that under case-cohort designs it is common to collect second-phase covariate data for more than a single outcome, since the case-cohort design naturally enables the analysis of multiple outcomes using a single subcohort selection. This is also the case in our example cohort discussed in Sections1 and 6.3, where, in addition to cases of all-cause mortality, genotype data has been collected also on cases of nonfatal incident cardiovascular disease events. To keep the example simple, in the data analysis, we consider only the case-cohort set for all-cause mortality. However, it is in principle straightforward to accommodate multiple outcome types in the likelihood expressions discussed below by using competing risks type notation where Ti denotes the observed time of the first incident disease event, death, or censoring, with Ei=0 indicating censoring and Ei∈{1,2,…,K} the different types of outcome events for which the second-phase covariate data is collected. Since utilizing the likelihood expressions does not necessitate estimation of K different hazard functions, the endpoint definitions may be pooled as seen suitable.


#### 5.2.1. Full Likelihood

The full likelihood expression (2.1) is now conditioned on the selection rule T𝒞≥b𝒞≡Ti≥bi for all  i∈𝒞, that is, event-free survival until the age at the cohort baseline. The likelihood expression becomes   (5.4)P(RC,TC,EC,YC,ZO∣TC≥bC,XC,θ)∝(α,β,γ)∏i∈OP(Ti,Ei∣Xi,Yi,Zi,α)P(Yi∣Xi,Zi,β)P(Zi∣Xi,γ)∫yi∫ziP(Ti≥bi∣Xi,yi,zi,α)P(dyi∣Xi,zi,β)P(dzi∣Xi,γ)×∏i∈C∖O∫ziP(Ti,Ei∣Xi,Yi,zi,α)P(Yi∣Xi,zi,β)P(dzi∣Xi,γ)∫yi∫ziP(Ti≥bi∣Xi,yi,zi,α)P(dyi∣Xi,zi,β)P(dzi∣Xi,γ).


#### 5.2.2. Inverse Probability Weighting

The weighted pseudolikelihood approach does not readily take into account additional selection which occurs in studies of incident outcomes. This was also pointed out by Reilly et al. [15], who in their discussion note that their reweighting approach for case-control studies is not valid for incident or nested case-control studies, since the incident cases and healthy controls available cannot be reweighted to represent the general population. Although this is true for the basic estimating function (3.1), similarly as in the missing data case, a double weighting mechanism can be devised to account for the additional selection step as   (5.5)∑i∈Olog⁡P(Yi∣Xi,Zi,β)P(Ti≥bi∣Xi,Yi,Zi,α̂)P(Ri=1∣TC,EC,XC,YC).
Since this weighting requires estimation of the time-to-event model parameters α, resampling-based variance estimation would be required to obtain valid standard errors for the estimates of β to account for the uncertainty in estimation of the weights.


#### 5.2.3. Conditional Likelihood

As in Section4.2, we get   (5.6)P(QO∣RC,TC≥bC,QC∖O,XC,ZO,θ)=∏i∈OP(Qi∣Xi,Zi,θ)P(RC,TC≥bC∣QC∖O,XC,ZO,θ),
where the correction term becomes   (5.7)P(RC,TC≥bC∣QC∖O,XC,ZO,θ)=∫yi:i∈O∫ti∈[bi,∞):i∈O∑ei:i∈OP(RC∣TC,EC,YC,XC)×∏i∈O[(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)].
In the case-cohort design discussed in Section 4.3.1, with type I censoring, this further simplifies into   (5.8)P(RC,TC≥bC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)∝(α,β)∏i∈O∫yi∫ti∈[bi,ci]∑eiP(Ri=1∣ti,ei,yi,Xi)×P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)=∏i∈O∫yi[P(bi≤T̃i<ci∣Xi,yi,Zi,α)+π(bi)P(T̃i≥ci∣Xi,yi,Zi,α)]×P(Yi∈dyi∣Xi,Zi,β)=∏i∈O∫yi[P(T̃i≥bi∣Xi,yi,Zi,α)-(1-π(bi))P(T̃i≥ci∣Xi,yi,Zi,α)]×P(Yi∈dyi∣Xi,Zi,β).


## 5.1. Missing Second-Phase Covariate Data

Further issue to be considered in practical applications is possible missing covariate data within the set𝒪. For instance, if the covariates to be collected under the cohort sampling design are genotypes, after selection of subjects to be genotyped, it may turn out that the DNA amount or concentration is too low or the genotyping is otherwise unsuccessful. Let M𝒪≡{Mi:i∈𝒪} be a set of indicator variables indicating whether the measurement of Zi turned out to be unsuccessful after selection of subject i into the set 𝒪, and let ℳ≡{i:Mi=1}⊆𝒪. Now Zi is observed only on the set 𝒪∖ℳ. We consider the implications of this separately for all of the three approaches discussed above.


### 5.1.1. Full Likelihood

If the missingness can be assumed missing at random, that is,M𝒪⊥(Zℳ,θ)∣R𝒞,T𝒞,E𝒞,X𝒞,Y𝒞,Z𝒪∖ℳ, the missingness indicators can again be included in the proportionality constant, and the full likelihood expression becomes (5.1)P(RC,MO,TC,EC,YC,ZO∖M∣XC,θ)∝(α,β,γ)∏i∈O∖MP(Ti,Ei∣Xi,Yi,Zi,α)P(Yi∣Xi,Zi,β)P(Zi∣Xi,γ)×∏i∈C∖O∪M∫ziP(Ti,Ei∣Xi,Yi,zi,α)P(Yi∣Xi,zi,β)P(dzi∣Xi,γ).


### 5.1.2. Inverse Probability Weighting

If the missingness mechanism would be known (and missing at random), a weighted pseudolikelihood expression for the set𝒪∖ℳ could be obtained by multiplying the original selection probabilities by the probabilities of observing the value of Zi after selection   (5.2)∑i∈O∖Mlog⁡P(Yi∣Xi,Zi,β)P(Mi=0∣RC,TC,EC,XC,YC,ZO∖M)P(Ri=1∣TC,EC,XC,YC).
In practice, however, the missingness mechanism would have to be modeled to obtain estimates for these probabilities.


### 5.1.3. Conditional Likelihood

By partitioning the observed data intoW=(T𝒪∖ℳ,E𝒪∖ℳ,Y𝒪∖ℳ) and V=(R𝒞,M𝒪,T𝒞∖𝒪∪ℳ,E𝒞∖𝒪∪ℳ,Y𝒞∖𝒪∪ℳ,X𝒞,Z𝒪∖ℳ), we obtain a conditional likelihood (5.3)P(QO∖M∣RC,MO,QC∖O∪M,XC,ZO∖M,θ)∝θP(MO∣RC,QC,XC,ZO∖M,θ)P(MO∣RC,QC∖O∪M,XC,ZO∖M,θ)∏i∈O∖MP(Qi∣Xi,Zi,θ)P(RC∣QC∖O∪M,XC,ZO∖M,θ),
where P(R𝒞∣Q𝒞∖𝒪∪ℳ,X𝒞,Z𝒪∖ℳ,θ) is obtained from (4.5) by replacing the set 𝒪 with 𝒪∖ℳ. Unlike the second-phase sampling mechanism, the missingness mechanism is generally unknown. From (5.3), it can be seen that if we are willing to assume that Zℳ are either missing completely at random or that M𝒪 depends only on X𝒞, the terms involving the missingness indicators cancel out of the likelihood. On the other hand, if the missingness may depend on the response variables Q𝒞, the missingness mechanism would have to be modeled. However, in such a case, it may be easier to work with the partitioning W=(M𝒪,T𝒪,E𝒪,Y𝒪,Z𝒪∖ℳ) and V=(R𝒞,T𝒞∖𝒪,E𝒞∖𝒪,Y𝒞∖𝒪,X𝒞) and model the population distribution of Zi rather than trying to estimate parameters describing the missingness mechanism. This in fact corresponds to what was done by Saarela and Kulathinal [19] and was required there because haplotypes are only partially identifiable from unphased genotype data.


## 5.1.1. Full Likelihood

If the missingness can be assumed missing at random, that is,M𝒪⊥(Zℳ,θ)∣R𝒞,T𝒞,E𝒞,X𝒞,Y𝒞,Z𝒪∖ℳ, the missingness indicators can again be included in the proportionality constant, and the full likelihood expression becomes (5.1)P(RC,MO,TC,EC,YC,ZO∖M∣XC,θ)∝(α,β,γ)∏i∈O∖MP(Ti,Ei∣Xi,Yi,Zi,α)P(Yi∣Xi,Zi,β)P(Zi∣Xi,γ)×∏i∈C∖O∪M∫ziP(Ti,Ei∣Xi,Yi,zi,α)P(Yi∣Xi,zi,β)P(dzi∣Xi,γ).


## 5.1.2. Inverse Probability Weighting

If the missingness mechanism would be known (and missing at random), a weighted pseudolikelihood expression for the set𝒪∖ℳ could be obtained by multiplying the original selection probabilities by the probabilities of observing the value of Zi after selection   (5.2)∑i∈O∖Mlog⁡P(Yi∣Xi,Zi,β)P(Mi=0∣RC,TC,EC,XC,YC,ZO∖M)P(Ri=1∣TC,EC,XC,YC).
In practice, however, the missingness mechanism would have to be modeled to obtain estimates for these probabilities.


## 5.1.3. Conditional Likelihood

By partitioning the observed data intoW=(T𝒪∖ℳ,E𝒪∖ℳ,Y𝒪∖ℳ) and V=(R𝒞,M𝒪,T𝒞∖𝒪∪ℳ,E𝒞∖𝒪∪ℳ,Y𝒞∖𝒪∪ℳ,X𝒞,Z𝒪∖ℳ), we obtain a conditional likelihood (5.3)P(QO∖M∣RC,MO,QC∖O∪M,XC,ZO∖M,θ)∝θP(MO∣RC,QC,XC,ZO∖M,θ)P(MO∣RC,QC∖O∪M,XC,ZO∖M,θ)∏i∈O∖MP(Qi∣Xi,Zi,θ)P(RC∣QC∖O∪M,XC,ZO∖M,θ),
where P(R𝒞∣Q𝒞∖𝒪∪ℳ,X𝒞,Z𝒪∖ℳ,θ) is obtained from (4.5) by replacing the set 𝒪 with 𝒪∖ℳ. Unlike the second-phase sampling mechanism, the missingness mechanism is generally unknown. From (5.3), it can be seen that if we are willing to assume that Zℳ are either missing completely at random or that M𝒪 depends only on X𝒞, the terms involving the missingness indicators cancel out of the likelihood. On the other hand, if the missingness may depend on the response variables Q𝒞, the missingness mechanism would have to be modeled. However, in such a case, it may be easier to work with the partitioning W=(M𝒪,T𝒪,E𝒪,Y𝒪,Z𝒪∖ℳ) and V=(R𝒞,T𝒞∖𝒪,E𝒞∖𝒪,Y𝒞∖𝒪,X𝒞) and model the population distribution of Zi rather than trying to estimate parameters describing the missingness mechanism. This in fact corresponds to what was done by Saarela and Kulathinal [19] and was required there because haplotypes are only partially identifiable from unphased genotype data.


## 5.2. Incident Outcomes and Left Truncation

Previous discussion was specific to a primary mortality outcome using time on study as the main time scale. In this section, we discuss separately how the different methods can accommodate cohort sampling for incident nonfatal primary outcomes. In the analysis of secondary, non-time-to-event outcomes, the presence of left truncation due to exclusion of cases of prevalent disease presents an additional complication. If the parameters of the secondary outcome model correspond to the background population alive at the cohort baseline (rather than to the disease-free population), this additional selection factor requires further adjustment. If the primary outcome is a mortality endpoint, this is not an issue, since then there is no further selection due to prevalent conditions. In likelihood-based adjustment for left truncation, the main time scale of the analysis has to be chosen as age instead of time on study. In survival modeling, it is well known that conditioning on event-free survival until the age at study baseline corresponds to exclusion of the followup time before that (e.g., [37, page 580], and the references therein). However, this is no longer true in the case of missing covariate data (e.g., [8, pages 5997-5998]), or indeed in the analysis of secondary outcomes, as will be demonstrated below. The set notations 𝒞 and 𝒪 are here taken to refer to the disease free cohort and second-phase study group.It should also be noted that under case-cohort designs it is common to collect second-phase covariate data for more than a single outcome, since the case-cohort design naturally enables the analysis of multiple outcomes using a single subcohort selection. This is also the case in our example cohort discussed in Sections1 and 6.3, where, in addition to cases of all-cause mortality, genotype data has been collected also on cases of nonfatal incident cardiovascular disease events. To keep the example simple, in the data analysis, we consider only the case-cohort set for all-cause mortality. However, it is in principle straightforward to accommodate multiple outcome types in the likelihood expressions discussed below by using competing risks type notation where Ti denotes the observed time of the first incident disease event, death, or censoring, with Ei=0 indicating censoring and Ei∈{1,2,…,K} the different types of outcome events for which the second-phase covariate data is collected. Since utilizing the likelihood expressions does not necessitate estimation of K different hazard functions, the endpoint definitions may be pooled as seen suitable.


### 5.2.1. Full Likelihood

The full likelihood expression (2.1) is now conditioned on the selection rule T𝒞≥b𝒞≡Ti≥bi for all  i∈𝒞, that is, event-free survival until the age at the cohort baseline. The likelihood expression becomes   (5.4)P(RC,TC,EC,YC,ZO∣TC≥bC,XC,θ)∝(α,β,γ)∏i∈OP(Ti,Ei∣Xi,Yi,Zi,α)P(Yi∣Xi,Zi,β)P(Zi∣Xi,γ)∫yi∫ziP(Ti≥bi∣Xi,yi,zi,α)P(dyi∣Xi,zi,β)P(dzi∣Xi,γ)×∏i∈C∖O∫ziP(Ti,Ei∣Xi,Yi,zi,α)P(Yi∣Xi,zi,β)P(dzi∣Xi,γ)∫yi∫ziP(Ti≥bi∣Xi,yi,zi,α)P(dyi∣Xi,zi,β)P(dzi∣Xi,γ).


### 5.2.2. Inverse Probability Weighting

The weighted pseudolikelihood approach does not readily take into account additional selection which occurs in studies of incident outcomes. This was also pointed out by Reilly et al. [15], who in their discussion note that their reweighting approach for case-control studies is not valid for incident or nested case-control studies, since the incident cases and healthy controls available cannot be reweighted to represent the general population. Although this is true for the basic estimating function (3.1), similarly as in the missing data case, a double weighting mechanism can be devised to account for the additional selection step as   (5.5)∑i∈Olog⁡P(Yi∣Xi,Zi,β)P(Ti≥bi∣Xi,Yi,Zi,α̂)P(Ri=1∣TC,EC,XC,YC).
Since this weighting requires estimation of the time-to-event model parameters α, resampling-based variance estimation would be required to obtain valid standard errors for the estimates of β to account for the uncertainty in estimation of the weights.


### 5.2.3. Conditional Likelihood

As in Section4.2, we get   (5.6)P(QO∣RC,TC≥bC,QC∖O,XC,ZO,θ)=∏i∈OP(Qi∣Xi,Zi,θ)P(RC,TC≥bC∣QC∖O,XC,ZO,θ),
where the correction term becomes   (5.7)P(RC,TC≥bC∣QC∖O,XC,ZO,θ)=∫yi:i∈O∫ti∈[bi,∞):i∈O∑ei:i∈OP(RC∣TC,EC,YC,XC)×∏i∈O[(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)].
In the case-cohort design discussed in Section 4.3.1, with type I censoring, this further simplifies into   (5.8)P(RC,TC≥bC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)∝(α,β)∏i∈O∫yi∫ti∈[bi,ci]∑eiP(Ri=1∣ti,ei,yi,Xi)×P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)=∏i∈O∫yi[P(bi≤T̃i<ci∣Xi,yi,Zi,α)+π(bi)P(T̃i≥ci∣Xi,yi,Zi,α)]×P(Yi∈dyi∣Xi,Zi,β)=∏i∈O∫yi[P(T̃i≥bi∣Xi,yi,Zi,α)-(1-π(bi))P(T̃i≥ci∣Xi,yi,Zi,α)]×P(Yi∈dyi∣Xi,Zi,β).


## 5.2.1. Full Likelihood

The full likelihood expression (2.1) is now conditioned on the selection rule T𝒞≥b𝒞≡Ti≥bi for all  i∈𝒞, that is, event-free survival until the age at the cohort baseline. The likelihood expression becomes   (5.4)P(RC,TC,EC,YC,ZO∣TC≥bC,XC,θ)∝(α,β,γ)∏i∈OP(Ti,Ei∣Xi,Yi,Zi,α)P(Yi∣Xi,Zi,β)P(Zi∣Xi,γ)∫yi∫ziP(Ti≥bi∣Xi,yi,zi,α)P(dyi∣Xi,zi,β)P(dzi∣Xi,γ)×∏i∈C∖O∫ziP(Ti,Ei∣Xi,Yi,zi,α)P(Yi∣Xi,zi,β)P(dzi∣Xi,γ)∫yi∫ziP(Ti≥bi∣Xi,yi,zi,α)P(dyi∣Xi,zi,β)P(dzi∣Xi,γ).


## 5.2.2. Inverse Probability Weighting

The weighted pseudolikelihood approach does not readily take into account additional selection which occurs in studies of incident outcomes. This was also pointed out by Reilly et al. [15], who in their discussion note that their reweighting approach for case-control studies is not valid for incident or nested case-control studies, since the incident cases and healthy controls available cannot be reweighted to represent the general population. Although this is true for the basic estimating function (3.1), similarly as in the missing data case, a double weighting mechanism can be devised to account for the additional selection step as   (5.5)∑i∈Olog⁡P(Yi∣Xi,Zi,β)P(Ti≥bi∣Xi,Yi,Zi,α̂)P(Ri=1∣TC,EC,XC,YC).
Since this weighting requires estimation of the time-to-event model parameters α, resampling-based variance estimation would be required to obtain valid standard errors for the estimates of β to account for the uncertainty in estimation of the weights.


## 5.2.3. Conditional Likelihood

As in Section4.2, we get   (5.6)P(QO∣RC,TC≥bC,QC∖O,XC,ZO,θ)=∏i∈OP(Qi∣Xi,Zi,θ)P(RC,TC≥bC∣QC∖O,XC,ZO,θ),
where the correction term becomes   (5.7)P(RC,TC≥bC∣QC∖O,XC,ZO,θ)=∫yi:i∈O∫ti∈[bi,∞):i∈O∑ei:i∈OP(RC∣TC,EC,YC,XC)×∏i∈O[(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)].
In the case-cohort design discussed in Section 4.3.1, with type I censoring, this further simplifies into   (5.8)P(RC,TC≥bC∣TC∖O,EC∖O,YC∖O,XC,ZO,θ)∝(α,β)∏i∈O∫yi∫ti∈[bi,ci]∑eiP(Ri=1∣ti,ei,yi,Xi)×P(Ti∈dti,Ei=ei∣Xi,yi,Zi,α)P(Yi∈dyi∣Xi,Zi,β)=∏i∈O∫yi[P(bi≤T̃i<ci∣Xi,yi,Zi,α)+π(bi)P(T̃i≥ci∣Xi,yi,Zi,α)]×P(Yi∈dyi∣Xi,Zi,β)=∏i∈O∫yi[P(T̃i≥bi∣Xi,yi,Zi,α)-(1-π(bi))P(T̃i≥ci∣Xi,yi,Zi,α)]×P(Yi∈dyi∣Xi,Zi,β).


## 6. Illustrations

### 6.1. Simulation Study

In order to compare the efficiency of the alternative estimation methods discussed above, namely, full likelihood, conditional likelihood, and weighted pseudolikelihood, we supplemented the cohort data described in Section1 with a simulated covariate, following the real-data-based simulation approach described by Saarela et al. in [8, pages 5998–6000], in order to have the simulation setting resemble as closely as possible the real data analysis setting of Section 6.3. This Bayesian procedure corresponds to multiple imputation [21], where given the observed data (T𝒞,E𝒞,X𝒞,Y𝒞) and predetermined associations with the primary and secondary outcomes, an additional “missing” covariate Z𝒞 is simulated from the posterior predictive distribution P(Z𝒞∣T𝒞,E𝒞,X𝒞,Y𝒞). Case-cohort sampling was then simulated in the complete datasets (T𝒞,E𝒞,X𝒞,Y𝒞,Z𝒞) thus obtained. The objective of this approach, as compared to using completely simulated data, was to obtain simulation results directly relevant to the real study of interest. We included a set of N=5039 individuals with a BMI measurement (Yi) available and eligible for the case-cohort study. The primary (time-to-event) outcome (Ti,Ei) was taken to be all-cause mortality, with d=996. Given the primary and secondary outcome data and age at baseline (Xi≡bi) observed for the cohort, and the current parameter values, additional binary covariate values were drawn from the conditional distributions (6.1)P(Zi=zi∣Ti,Ei,Yi,Xi,α,β,γ)=P(Ti,Ei∣Xi,Yi,zi,α)P(Yi∣Xi,zi,β)P(Zi=zi∣Xi,γ)∑zi∈{0,1}P(Ti,Ei∣Xi,Yi,zi,α)P(Yi∣Xi,zi,β)P(Zi=zi∣Xi,γ),
where the model for the survival outcome, P(Ti∈dti,Ei=ei∣Xi,Yi,Zi,α)∝α[λi(ti)]1{ei=1}Si(ti), was specified as a proportional hazards model λi(u)≡λ0(u)exp{α1bi+α2Yi+α3Zi}, with time on study as the time scale and age at baseline as a covariate. Here, the the baseline hazard function λ0 was specified in terms of a piecewise constant function using 15 time bins of equal length over the followup period of seven years. Flat improper priors were used for all parameters (those not fixed to selected values), on log scale for the nonnegative parameters.Normal modelYi∣Xi,Zi,β~N(β0+β1bi+β2Zi,β32) was used for the secondary outcome, while the population distribution of the additional covariate was specified as P(Zi=zi∣Xi,γ)=γzi(1-γ)1-zi. Here, parameters α3, β2, and γ were fixed at selected values while the other parameters were allowed to be determined by the data and integrated out by drawing from their respective full conditional posterior distributions using Markov chain Monte Carlo sampling. Variables Yi and bi were centered but otherwise untransformed. 1000 values for each Zi, i∈𝒞, were obtained by running the MCMC sampler for 25000 rounds after a 10000-round burn-in with each set of fixed values of (α3,β2,γ) given in Table 1 and saving every 25th state of the chain. In each of the complete datasets obtained by combining the observed data and the simulated covariate values Zi, case-cohort selection was carried out by matching the age distribution of the subcohort to that of the cases by first fitting a logistic regression model with age as a covariate to the observed survival status, as detailed in Section 4.3.1. The predictive probabilities from the logistic model were scaled to give an expected subcohort size of 1000, and the subcohort was selected using the obtained probabilities in Bernoulli sampling. This procedure gave an expected case-cohort set size of E(n∣E𝒞,X𝒞)=1761.Table 1
Maximum likelihood estimates of the parametersα of the all-cause mortality model λi(t)=(α4/α0)(t/α0)α4-1exp{α1bi+α2Yi+α3Zi}, β of the BMI model Yi∣Xi,Zi,β~N(β0+β1bi+β2Zi,β32), and γ, the population frequency of the simulated covariate (true value 0.25). Parameters of main interest are α2 and β2, which characterize the association between the simulated covariate and primary (all-cause mortality) and secondary (BMI) outcomes, respectively. The values are means and standard deviations over 1000 replications. SE stands for the standard errors of the maximum likelihood estimates.
α3β2Lik.α̂0α̂1α̂2α̂3SD(α̂3)SE(α̂3)α̂4β̂0β̂1β̂2SD(β̂2)SE(β̂2)β̂3γ̂0.00.0Full23.720.09-0.01-0.010.100.101.330.00-0.080.010.220.224.100.25Cond.23.740.09-0.010.000.100.101.330.00-0.080.000.230.234.17—Weighted———————-0.01-0.080.030.300.304.10—0.50.0Full26.300.09-0.010.490.090.091.330.00-0.080.010.220.214.100.25Cond.26.330.09-0.010.500.090.091.330.00-0.080.000.220.224.17—Weighted———————-0.01-0.080.030.300.304.10—0.00.5Full23.730.09-0.01-0.010.100.101.33-0.13-0.080.530.250.234.090.25Cond.23.750.09-0.010.000.100.101.33-0.13-0.080.520.240.234.17—Weighted———————-0.13-0.080.530.310.314.09—0.50.5Full26.320.09-0.020.490.090.091.33-0.14-0.080.550.230.224.090.25Cond.26.340.09-0.010.500.090.091.33-0.14-0.080.540.220.224.17—Weighted———————-0.14-0.080.540.310.314.09—1.00.0Full29.620.09-0.010.990.090.091.350.00-0.080.020.210.214.100.25Cond.29.640.09-0.011.000.090.091.350.00-0.080.010.220.214.17—Weighted———————-0.01-0.080.030.290.294.10—0.01.0Full23.770.09-0.010.000.100.101.33-0.28-0.081.110.260.234.070.25Cond.23.770.09-0.010.000.100.101.33-0.27-0.081.060.240.234.15—Weighted———————-0.26-0.081.050.320.324.07—1.01.0Full29.670.09-0.031.000.090.091.35-0.29-0.081.140.240.224.070.25Cond.29.650.09-0.021.000.090.091.35-0.28-0.081.080.220.214.14—Weighted———————-0.27-0.081.070.310.314.07—In fitting models to the datasets so obtained, we used the same model specifications as above, with the exception of the proportional hazards model, where we fitted a “misspecified” Weibull modelλi(u)≡(α4/α0)(u/α0)α4-1exp{α1bi+α2Yi+α3Zi}. Even though this is a different model than the one specified for the simulation step, we did not simulate the time-to-event outcome data, which was taken from the example cohort; a comparison between fitted piecewise constant and Weibull hazards (results not shown) indicated that the Weibull model is adequate for modeling these data, accounting for the increasing baseline mortality rate due to ageing of the cohort. Full likelihood function (2.1), conditional likelihood function (4.3), and weighted pseudolikelihood function (3.1) were maximized with respect to the parameters (α,β,γ), (α,β), and β, respectively, substituting in the above parametric models. The numerical optimization was carried out using the optim function of the R statistical software, applying the BFGS optimization algorithm [38]. The integrals over yi in (4.7) were evaluated using numerical (quadrature) integration, using the gsl_integration_qagi function of GNU scientific library [39]. The standard errors of the maximum likelihood estimators were evaluated by inverting the numerically differentiated Hessian matrix at the maximum likelihood point. The standard errors of the inverse-probability-weighted estimates were evaluated using the robust variance formula of Appendix A with the observed information (A.4) and observed score (A.5). The results from the 1000 replications are presented in Table 1. Here, the parameters (α3,β2,γ), corresponding to the simulated covariate, are of main interest; the other parameters reflect the observed data, and thus the Monte Carlo variances for these are not relevant. The standard errors for the corresponding parameters are given in the real-data example of Section 6.3. The results indicate that conditional likelihood estimation gives similar efficiency compared to the corresponding full likelihood estimates irrespective of whether the covariate of interest was associated with the survival outcome. In contrast, both likelihood-based approaches gave better efficiency compared to inverse probability weighting. This is expected, as the weighted pseudolikelihood gives smaller weights to cases, whereas the conditional likelihood weights cases and noncases differentially only to the extent they actually differ. Both the inverse observed information and robust variance estimates agreed well with the Monte Carlo variances.


### 6.2. Multimodality under Full Likelihood When the Sampling Fraction Is Small

Let now the observed data be only(R𝒞,Y𝒞,Z𝒪), and let the sampling mechanism be simple random sampling without replacement; P(R𝒞)=1/(Nn), where the sample size n is fixed. The specific aim of the following example is to demonstrate the vulnerability of observed data likelihoods, integrated over the missing covariate data on the set 𝒞∖𝒪, to misspecification of the response model. The full likelihood expression now becomes P(R𝒞,Y𝒞,Z𝒪∣θ)∝(β,γ)∏i∈𝒪[P(Yi∣Zi,β)P(Zi∣γ)]∏i∈𝒞∖𝒪∫ziP(Yi∣zi,β)P(Zi∈dzi∣γ). Due to the simple random sampling, the corresponding conditional likelihoods simplify into P(Y𝒪∣R𝒞,Y𝒞∖𝒪,Z𝒪,θ)=∏i∈𝒪P(Yi∣Zi,β) and P(Y𝒪,Z𝒪∣R𝒞,Y𝒞∖𝒪,θ)=∏i∈𝒪[P(Yi∣Zi,β)P(Zi∣γ)].WithN=1000, we simulated covariate values from Zi~Bernoulli(0.2) and two different sets of response values (both independently of Zi) from Yi~N(4,4) and Yi~Gamma(4,1) (same mean and variance, but the latter distribution is skewed to the right). We fit two alternative models to these data, in both models Zi∣γ~Bernoulli(γ), with the models for the response specified as Yi∣Zi,β~N(β0+β1Zi,β22) and Yi∣Zi,β~NPQM(β0+β1Zi,β2,β3,β4), where β are the collections of all model parameters. The latter model is the normal-polynomial quantile mixture distribution proposed by Karvanen ([20, pages 950–953]; see also R package Lmoments), which we apply here for regression modeling. The parametrization here can be expressed in terms of the first four L-moments as λ1=E(Yi∣Zi,β)=β0+β1Zi, λ2=β2/π (L-scale), λ3=β3 (L-skewness), and λ4=β4 (L-kurtosis). This distribution is suitable for regression modeling as the first L-moment is the mean of the distribution, while the two additional shape parameters allow for more flexible modeling of the residual distribution. The NPQM distribution includes normal distribution as a special case when λ3=0 and λ4/λ2=0.1226. The full and conditional log-likelihood expressions under these model specifications were maximized with respect to β and γ as described in the previous section. Initial values for the algorithm were set at the sample moments calculated from the marginal distribution of the response variable, with the regression coefficient set to β1=0 and γ=0.5.Table2 shows the mean maximum likelihood estimates over 1000 replications for all 8 combinations of data generating model, fitted model, and estimation method. With the response simulated from normal distribution, both models estimated using either full or conditional likelihood give the expected results (Figures 1 and 2). However, with the response simulated from gamma distribution and the sampling fraction n/N small, the misspecified normal model estimated using full likelihood indicates a spurious association between the response and the covariate. This is because the missing data act as extra parameters, allowing the optimization procedure to obtain a better fit to the skewed data. The corresponding sampling distributions for β1 shown in Figure 3 have become bimodal. It should be noted that the initial values given to the optimization algorithm corresponded always to the “correct solution” of no covariate effect, thus enabling the algorithm to find the correct mode. However, bimodality in the sampling distribution does not necessarily indicate that the likelihood functions given a single realized dataset would be bimodal.Table 2
Maximum likelihood estimates of parametersβ of the outcome models Yi∣Zi,β~N(β0+β1Zi,β22) and Yi∣Zi,β~NPQM(β0+β1Zi,β2,β3,β4) and the population frequency γ of the binary covariate Zi (means over 1000 replications). Outcome variable was simulated from Yi~N(4,4) and Yi~Gamma(4,1).
Normal datanNormal modelNPQM modelβ̂0β̂1SE(β̂1)β̂2γ̂β̂0β̂1SE(β̂1)β̂2β̂3SE(β̂3)β̂4γ̂Full lik.10004.000.000.162.000.204.000.000.162.000.000.010.140.205004.000.010.221.990.204.000.010.221.990.000.010.140.202004.000.000.351.990.204.000.000.351.990.000.010.140.201004.000.010.501.990.204.000.010.501.990.000.020.140.20Cond. lik.5004.000.010.221.990.204.000.010.221.990.000.020.140.202004.010.000.351.990.204.010.000.351.990.000.030.130.201004.000.010.501.980.204.000.010.501.980.000.050.130.20Gamma datanNormal modelNPQM modelβ̂0β̂1SE(β̂1)β̂2γ̂β̂0β̂1SE(β̂1)β̂2β̂3SE(β̂3)β̂4γ̂Full lik.10004.000.010.161.990.203.990.000.132.010.150.010.140.205004.000.010.221.990.203.990.010.202.010.150.010.140.202003.970.150.381.980.203.980.070.352.000.150.020.140.201003.860.910.381.880.193.950.180.321.960.160.020.140.22Cond. lik.5004.000.000.221.990.204.000.000.192.010.150.020.140.202004.00-0.010.351.980.204.010.000.292.000.150.030.140.201004.00-0.010.501.970.204.010.010.412.000.150.050.130.20Figure 1
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; normal model fitted to normal data (Yi~N(4,4)).Figure 2
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; NPQM model fitted to normal data (Yi~N(4,4)).Figure 3
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; normal model fitted to gamma data (Yi~Gamma(4,1)).Using the NPQM model which allows a skewed residual distribution does not correct the situation either when combined with full likelihood estimation (Figure4). This is not unexpected since adding more parameters to an already overparameterized situation is not necessarily helpful in solving identification problems. In contrast, conditional likelihood continues to give reasonable results in all cases, with the skewness of the data correctly reflected by the parameter β3 of the NPQM model. Some degree of misspecification of the response model is required for the multimodality to appear since when we simulated the response variable from NPQM distribution, the full likelihood fit of the NPQM model indicated no problems, while the full likelihood fit of the normal model showed the same problems as with gamma data (Figures 5 and 6). We also repeated all the simulation alternatives with β1=1 (Yi~Zi+N(4,4) and Yi~Zi+Gamma(4,1)), with the conclusions essentially unchanged (Table 3 and Figures 7, 8, 9, and 10).Table 3
Maximum likelihood estimates of parametersβ of the outcome models Yi∣Zi,β~N(β0+β1Zi,β22) and Yi∣Zi,β~NPQM(β0+β1Zi,β2,β3,β4) and the population frequency γ of the binary covariate Zi (means over 1000 replications). Outcome variable was simulated from Yi∣Zi~Zi+N(4,4) and Yi∣Zi~Zi+Gamma(4,1).
Normal datanNormal modelNPQM modelβ̂0β̂1SE(β̂1)β̂2γ̂β̂0β̂1SE(β̂1)β̂2β̂3SE(β̂3)β̂4γ̂Full lik.10004.001.000.162.000.204.001.000.162.000.000.010.140.205004.001.010.221.990.204.001.010.221.990.000.020.140.202004.001.000.341.990.204.001.000.341.990.000.020.140.201004.001.000.471.990.204.001.010.481.990.000.020.140.20Cond. lik.5004.001.010.221.990.204.001.010.221.990.000.020.140.202004.011.000.351.990.204.011.000.351.990.000.030.130.201004.001.010.501.980.204.001.010.511.980.000.050.130.20Gamma datanNormal modelNPQM modelβ̂0β̂1SE(β̂1)β̂2γ̂β̂0β̂1SE(β̂1)β̂2β̂3SE(β̂3)β̂4γ̂Full lik.10004.001.010.161.990.203.991.000.132.010.150.010.140.205003.971.150.241.980.203.971.080.182.000.150.020.140.202003.832.010.361.860.193.921.280.231.970.160.020.140.211003.722.950.261.700.163.871.440.261.940.160.020.140.23Cond. lik.5004.001.000.221.990.204.001.000.192.010.150.020.140.202004.000.990.351.980.204.011.000.292.000.150.030.140.201004.000.990.501.970.204.011.010.412.000.150.050.130.20Figure 4
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; NPQM model fitted to gamma data (Yi~Gamma(4,1)).Figure 5
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; normal model fitted to NPQM data (Yi~NPQM(4,2,0.15,0.14)).Figure 6
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; NPQM model fitted to NPQM data (Yi~NPQM(4,2,0.15,0.14)).Figure 7
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; normal model fitted to normal data (Yi~Zi+N(4,4)).Figure 8
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; NPQM model fitted to normal data (Yi~Zi+N(4,4)).Figure 9
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; normal model fitted to gamma data (Yi~Zi+Gamma(4,1)).Figure 10
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; NPQM model fitted to gamma data (Yi~Zi+Gamma(4,1)).


### 6.3. An Example with Real Data

The case-cohort set for all-cause mortality in the example cohort (N=5039) is of size n=1816, the union of a subcohort of size 1068 and 996 deaths due to any cause (since the subcohort has been selected independently of the case status, it includes a number of cases). The case-cohort design applied in this cohort was described in Section 4.3.1. Due to sample unavailability or unsuccessful genotyping, the genotype for the lactase persistence SNP rs4988235 was unavailable for 156 individuals, giving |𝒪∖ℳ|=1660. The model for the data described in Section 1 was specified as follows: Hardy-Weinberg model Zi∣γ~Binom(2,γ) was used for the number of lactase persistence alleles. The model for BMI, which is the response variable of interest, given age at baseline Xi≡bi and the genotype Zi, was taken to be Yi∣bi,Zi,β~NPQM(β0+β1bi+β21{Zi=0},β3,β4,β5). The model for the survival outcome was taken to be P(Ti∈dti,Ei=ei∣Yi,Zi,α)∝α[λi(ti)]1{ei=1}Si(ti), where Si(ti)≡exp{∫0tiλi(u)du}, and using again the Weibull form for the hazard function λi(u)≡(α4/α0)(u/α0)α4-1exp{α1bi+α2Yi+α31{Zi=0}}. In both regression models, the lactase persistence allele noncarriers are compared to the hetero- and homozygote carriers of the allele. Under the case-cohort design where the subcohort sampling probabilities depend on bi, the conditional likelihood expression takes the form   (6.2)P(TO∖M,EO∖M,YO∖M∣RC,MO,TC∖O∪M,EC∖O∪M,YC∖O∪M,XC,ZO∖M,θ)∝(α,β)∏i∈O∖MP(Ti,Ei∣Yi,Zi,α)P(Yi∣bi,Zi,β)∫yi[1-(1-π(bi))Si(ci)]P(dyi∣bi,Zi,β),
or   (6.3)P(MO,TO,EO,YO,ZO∖M∣RC,TC∖O,EC∖O,YC∖O,XC,θ)∝(α,β,γ)∏i∈O∖MP(Ti,Ei∣Yi,Zi,α)P(Yi∣bi,Zi,β)P(Zi=zi∣γ)∑zi∫yi[1-(1-π(bi))Si(ci)]P(dyi∣bi,zi,β)P(Zi=zi∣γ)×∏i∈M∑ziP(Ti,Ei∣Yi,zi,α)P(Yi∣bi,zi,β)P(Zi=zi∣γ)∑zi∫yi[1-(1-π(bi))Si(ci)]P(dyi∣bi,zi,β)P(Zi=zi∣γ),
depending on whether we condition upon the observed genotype data Z𝒪∖ℳ, or whether this is modeled as part of the likelihood (see Section 5.1.3). The difference between these two approaches is that the former requires more assumptions on the missingness mechanism. Therefore, we compare here the two approaches to see whether there is any observable difference between the results. In addition, expression (6.3) enables estimation of the population allele frequency γ. However, it should be noted that while the summation over the single SNP genotype variable Zi is not a computational problem in the present example, this is not necessarily the case generally, with multiple continuous covariates in the model. Therefore, the potential advantage of expression such as (6.2) is that it avoids specification of the covariate distribution and the integration over the missing covariate values.A remaining issue to be considered is the numerical evaluation of the integrals overyi in the correction terms of (6.2) and (6.3). It should be noted that these have to be evaluated at each parameter value tried at the hill-climbing numerical optimization. Although numerical (quadrature) integration would be feasible in one dimension, here we propose evaluation of such integrals using Monte Carlo integration by repeated sampling from the distribution P(Yi∈dyi∣bi,Zi,β). This has an added advantage that while the numerical evaluation of the density of the NPQM model is computationally expensive, random variates can be easily simulated from this distribution since it is defined through the quantile function. This proceeds by drawing random variates from [0,1]-uniform distribution and applying formula (6.2) of Karvanen [20]. With random variates yi(j)~NPQM(β0+β1bi+β21{Zi=0},β3,β4,β5), j=1,…,k, drawn for each i∈𝒪∖ℳ given the current parameter values, the integrals in (6.2) are then approximated by the means (1/k)∑j=1k[1-(1-π(bi))P(Ti>ci∣yi(j),Zi,α)].Maximum likelihood estimates obtained by numerical maximization of the expressions (6.2) and (6.3) with respect to α, β, and in the latter case also γ, are presented in Table 4 (α̂ and γ̂) and Table 5 (β̂). For comparison, we also maximised a full likelihood of the form (5.1) with respect to (α,β,γ) and a weighted pseudolikelihood of the form (5.2) with respect to β. In the latter case, we accounted for the missingness within the case-cohort set by fitting a logistic regression model logit{P(Mi=0∣Ei,bi,Yi,η)}=η0+η1Ei+η2bi+η3Yi in the set 𝒪 and calculated the adjusted weights as inverses of P(Mi=0∣Ei,bi,Yi,η̂)P(Ri=1∣Ei,bi). Due to the extra estimation step, we used bootstrap with 5000 replications to obtain standard errors for the weighted pseudolikelihood estimates. Both conditional likelihood expressions gave very similar results, suggesting that the missing data within the case-cohort set is not a major issue in the present case. The comparison between different number of Monte Carlo replicates in numerical evaluation of the conditional likelihood correction term suggests that the estimates do no longer appreciably change when k is increased from 1000. However, the change in the estimates when increasing k from 100 to 1000 is already well within the standard errors.Table 4
Maximum likelihood estimates (standard errors) of the parametersα of the all-cause mortality model λi(t)=(α4/α0)(t/α0)α4-1exp{α1bi+α2Yi+α31{Zi=0}}, and γ, the population frequency of the lactase persistence allele.
Likelihoodkα̂0logα̂0α̂1α̂2α̂3α̂4logα̂4γ̂Conditional10023.883.17 (0.21)0.12 (0.01)0.00 (0.01)-0.01 (0.12)1.340.29 (0.03)—(6.2)100023.523.16 (0.21)0.12 (0.01)0.00 (0.01)-0.01 (0.12)1.340.29 (0.03)—500023.453.15 (0.21)0.12 (0.01)0.00 (0.01)-0.01 (0.12)1.340.29 (0.03)—Conditional10020.433.02 (0.20)0.12 (0.01)-0.01 (0.01)-0.01 (0.11)1.330.29 (0.03)0.59 (0.01)(6.3)100020.343.01 (0.20)0.12 (0.01)-0.01 (0.01)-0.01 (0.11)1.330.29 (0.03)0.59 (0.01)500020.443.02 (0.20)0.12 (0.01)-0.01 (0.01)-0.01 (0.11)1.330.29 (0.03)0.59 (0.01)Full—18.402.91 (0.17)0.09 (0.01)-0.01 (0.01)-0.07 (0.11)1.330.28 (0.03)0.59 (0.01)Table 5
Maximum likelihood estimates (standard errors) of the parametersβ of the BMI model Yi∣bi,Zi,β~NPQM(β0+β1bi+β21{Zi=0},β3,β4,β5).
Likelihoodkβ̂0β̂1β̂2β̂3log(β̂3/π)β̂4β̂5logβ̂5Conditional10026.91 (0.14)-0.07 (0.02)-0.64 (0.27)4.270.88 (0.02)0.14 (0.03)0.37-1.00 (0.05)(6.2)100026.92 (0.14)-0.07 (0.02)-0.64 (0.27)4.270.88 (0.02)0.14 (0.03)0.37-1.00  (0.05)500026.92 (0.14)-0.07 (0.02)-0.64 (0.27)4.270.88 (0.02)0.14 (0.03)0.37-1.00 (0.05)Conditional10026.93 (0.13)-0.07 (0.02)-0.63 (0.26)4.250.87 (0.02)0.14 (0.03)0.37-1.00 (0.05)(6.3)100026.93 (0.13)-0.07 (0.02)-0.63 (0.26)4.250.87 (0.02)0.14 (0.03)0.37-1.00 (0.05)500026.93 (0.13)-0.07 (0.02)-0.63 (0.26)4.250.87 (0.02)0.14 (0.03)0.37-1.00 (0.05)Full—26.94 (0.07)-0.07 (0.01)-0.60 (0.25)4.060.83 (0.01)0.13 (0.01)0.35-1.06 (0.03)Weighted—26.94 (0.15)-0.06 (0.03)-0.68 (0.33)4.300.88 (0.04)0.17 (0.05)0.41-0.89 (0.11)Full likelihood and weighted pseudolikelihood estimates agreed well with the conditional likelihood ones, although the latter had higher standard errors, as was the case also in the simulations. As noted by Kettunen et al. [10], the absence of the lactase persistence allele shows association with lower body mass index. The residuals from the model for BMI are significantly skewed to the right, as indicated by the estimates of β4.


## 6.1. Simulation Study

In order to compare the efficiency of the alternative estimation methods discussed above, namely, full likelihood, conditional likelihood, and weighted pseudolikelihood, we supplemented the cohort data described in Section1 with a simulated covariate, following the real-data-based simulation approach described by Saarela et al. in [8, pages 5998–6000], in order to have the simulation setting resemble as closely as possible the real data analysis setting of Section 6.3. This Bayesian procedure corresponds to multiple imputation [21], where given the observed data (T𝒞,E𝒞,X𝒞,Y𝒞) and predetermined associations with the primary and secondary outcomes, an additional “missing” covariate Z𝒞 is simulated from the posterior predictive distribution P(Z𝒞∣T𝒞,E𝒞,X𝒞,Y𝒞). Case-cohort sampling was then simulated in the complete datasets (T𝒞,E𝒞,X𝒞,Y𝒞,Z𝒞) thus obtained. The objective of this approach, as compared to using completely simulated data, was to obtain simulation results directly relevant to the real study of interest. We included a set of N=5039 individuals with a BMI measurement (Yi) available and eligible for the case-cohort study. The primary (time-to-event) outcome (Ti,Ei) was taken to be all-cause mortality, with d=996. Given the primary and secondary outcome data and age at baseline (Xi≡bi) observed for the cohort, and the current parameter values, additional binary covariate values were drawn from the conditional distributions (6.1)P(Zi=zi∣Ti,Ei,Yi,Xi,α,β,γ)=P(Ti,Ei∣Xi,Yi,zi,α)P(Yi∣Xi,zi,β)P(Zi=zi∣Xi,γ)∑zi∈{0,1}P(Ti,Ei∣Xi,Yi,zi,α)P(Yi∣Xi,zi,β)P(Zi=zi∣Xi,γ),
where the model for the survival outcome, P(Ti∈dti,Ei=ei∣Xi,Yi,Zi,α)∝α[λi(ti)]1{ei=1}Si(ti), was specified as a proportional hazards model λi(u)≡λ0(u)exp{α1bi+α2Yi+α3Zi}, with time on study as the time scale and age at baseline as a covariate. Here, the the baseline hazard function λ0 was specified in terms of a piecewise constant function using 15 time bins of equal length over the followup period of seven years. Flat improper priors were used for all parameters (those not fixed to selected values), on log scale for the nonnegative parameters.Normal modelYi∣Xi,Zi,β~N(β0+β1bi+β2Zi,β32) was used for the secondary outcome, while the population distribution of the additional covariate was specified as P(Zi=zi∣Xi,γ)=γzi(1-γ)1-zi. Here, parameters α3, β2, and γ were fixed at selected values while the other parameters were allowed to be determined by the data and integrated out by drawing from their respective full conditional posterior distributions using Markov chain Monte Carlo sampling. Variables Yi and bi were centered but otherwise untransformed. 1000 values for each Zi, i∈𝒞, were obtained by running the MCMC sampler for 25000 rounds after a 10000-round burn-in with each set of fixed values of (α3,β2,γ) given in Table 1 and saving every 25th state of the chain. In each of the complete datasets obtained by combining the observed data and the simulated covariate values Zi, case-cohort selection was carried out by matching the age distribution of the subcohort to that of the cases by first fitting a logistic regression model with age as a covariate to the observed survival status, as detailed in Section 4.3.1. The predictive probabilities from the logistic model were scaled to give an expected subcohort size of 1000, and the subcohort was selected using the obtained probabilities in Bernoulli sampling. This procedure gave an expected case-cohort set size of E(n∣E𝒞,X𝒞)=1761.Table 1
Maximum likelihood estimates of the parametersα of the all-cause mortality model λi(t)=(α4/α0)(t/α0)α4-1exp{α1bi+α2Yi+α3Zi}, β of the BMI model Yi∣Xi,Zi,β~N(β0+β1bi+β2Zi,β32), and γ, the population frequency of the simulated covariate (true value 0.25). Parameters of main interest are α2 and β2, which characterize the association between the simulated covariate and primary (all-cause mortality) and secondary (BMI) outcomes, respectively. The values are means and standard deviations over 1000 replications. SE stands for the standard errors of the maximum likelihood estimates.
α3β2Lik.α̂0α̂1α̂2α̂3SD(α̂3)SE(α̂3)α̂4β̂0β̂1β̂2SD(β̂2)SE(β̂2)β̂3γ̂0.00.0Full23.720.09-0.01-0.010.100.101.330.00-0.080.010.220.224.100.25Cond.23.740.09-0.010.000.100.101.330.00-0.080.000.230.234.17—Weighted———————-0.01-0.080.030.300.304.10—0.50.0Full26.300.09-0.010.490.090.091.330.00-0.080.010.220.214.100.25Cond.26.330.09-0.010.500.090.091.330.00-0.080.000.220.224.17—Weighted———————-0.01-0.080.030.300.304.10—0.00.5Full23.730.09-0.01-0.010.100.101.33-0.13-0.080.530.250.234.090.25Cond.23.750.09-0.010.000.100.101.33-0.13-0.080.520.240.234.17—Weighted———————-0.13-0.080.530.310.314.09—0.50.5Full26.320.09-0.020.490.090.091.33-0.14-0.080.550.230.224.090.25Cond.26.340.09-0.010.500.090.091.33-0.14-0.080.540.220.224.17—Weighted———————-0.14-0.080.540.310.314.09—1.00.0Full29.620.09-0.010.990.090.091.350.00-0.080.020.210.214.100.25Cond.29.640.09-0.011.000.090.091.350.00-0.080.010.220.214.17—Weighted———————-0.01-0.080.030.290.294.10—0.01.0Full23.770.09-0.010.000.100.101.33-0.28-0.081.110.260.234.070.25Cond.23.770.09-0.010.000.100.101.33-0.27-0.081.060.240.234.15—Weighted———————-0.26-0.081.050.320.324.07—1.01.0Full29.670.09-0.031.000.090.091.35-0.29-0.081.140.240.224.070.25Cond.29.650.09-0.021.000.090.091.35-0.28-0.081.080.220.214.14—Weighted———————-0.27-0.081.070.310.314.07—In fitting models to the datasets so obtained, we used the same model specifications as above, with the exception of the proportional hazards model, where we fitted a “misspecified” Weibull modelλi(u)≡(α4/α0)(u/α0)α4-1exp{α1bi+α2Yi+α3Zi}. Even though this is a different model than the one specified for the simulation step, we did not simulate the time-to-event outcome data, which was taken from the example cohort; a comparison between fitted piecewise constant and Weibull hazards (results not shown) indicated that the Weibull model is adequate for modeling these data, accounting for the increasing baseline mortality rate due to ageing of the cohort. Full likelihood function (2.1), conditional likelihood function (4.3), and weighted pseudolikelihood function (3.1) were maximized with respect to the parameters (α,β,γ), (α,β), and β, respectively, substituting in the above parametric models. The numerical optimization was carried out using the optim function of the R statistical software, applying the BFGS optimization algorithm [38]. The integrals over yi in (4.7) were evaluated using numerical (quadrature) integration, using the gsl_integration_qagi function of GNU scientific library [39]. The standard errors of the maximum likelihood estimators were evaluated by inverting the numerically differentiated Hessian matrix at the maximum likelihood point. The standard errors of the inverse-probability-weighted estimates were evaluated using the robust variance formula of Appendix A with the observed information (A.4) and observed score (A.5). The results from the 1000 replications are presented in Table 1. Here, the parameters (α3,β2,γ), corresponding to the simulated covariate, are of main interest; the other parameters reflect the observed data, and thus the Monte Carlo variances for these are not relevant. The standard errors for the corresponding parameters are given in the real-data example of Section 6.3. The results indicate that conditional likelihood estimation gives similar efficiency compared to the corresponding full likelihood estimates irrespective of whether the covariate of interest was associated with the survival outcome. In contrast, both likelihood-based approaches gave better efficiency compared to inverse probability weighting. This is expected, as the weighted pseudolikelihood gives smaller weights to cases, whereas the conditional likelihood weights cases and noncases differentially only to the extent they actually differ. Both the inverse observed information and robust variance estimates agreed well with the Monte Carlo variances.


## 6.2. Multimodality under Full Likelihood When the Sampling Fraction Is Small

Let now the observed data be only(R𝒞,Y𝒞,Z𝒪), and let the sampling mechanism be simple random sampling without replacement; P(R𝒞)=1/(Nn), where the sample size n is fixed. The specific aim of the following example is to demonstrate the vulnerability of observed data likelihoods, integrated over the missing covariate data on the set 𝒞∖𝒪, to misspecification of the response model. The full likelihood expression now becomes P(R𝒞,Y𝒞,Z𝒪∣θ)∝(β,γ)∏i∈𝒪[P(Yi∣Zi,β)P(Zi∣γ)]∏i∈𝒞∖𝒪∫ziP(Yi∣zi,β)P(Zi∈dzi∣γ). Due to the simple random sampling, the corresponding conditional likelihoods simplify into P(Y𝒪∣R𝒞,Y𝒞∖𝒪,Z𝒪,θ)=∏i∈𝒪P(Yi∣Zi,β) and P(Y𝒪,Z𝒪∣R𝒞,Y𝒞∖𝒪,θ)=∏i∈𝒪[P(Yi∣Zi,β)P(Zi∣γ)].WithN=1000, we simulated covariate values from Zi~Bernoulli(0.2) and two different sets of response values (both independently of Zi) from Yi~N(4,4) and Yi~Gamma(4,1) (same mean and variance, but the latter distribution is skewed to the right). We fit two alternative models to these data, in both models Zi∣γ~Bernoulli(γ), with the models for the response specified as Yi∣Zi,β~N(β0+β1Zi,β22) and Yi∣Zi,β~NPQM(β0+β1Zi,β2,β3,β4), where β are the collections of all model parameters. The latter model is the normal-polynomial quantile mixture distribution proposed by Karvanen ([20, pages 950–953]; see also R package Lmoments), which we apply here for regression modeling. The parametrization here can be expressed in terms of the first four L-moments as λ1=E(Yi∣Zi,β)=β0+β1Zi, λ2=β2/π (L-scale), λ3=β3 (L-skewness), and λ4=β4 (L-kurtosis). This distribution is suitable for regression modeling as the first L-moment is the mean of the distribution, while the two additional shape parameters allow for more flexible modeling of the residual distribution. The NPQM distribution includes normal distribution as a special case when λ3=0 and λ4/λ2=0.1226. The full and conditional log-likelihood expressions under these model specifications were maximized with respect to β and γ as described in the previous section. Initial values for the algorithm were set at the sample moments calculated from the marginal distribution of the response variable, with the regression coefficient set to β1=0 and γ=0.5.Table2 shows the mean maximum likelihood estimates over 1000 replications for all 8 combinations of data generating model, fitted model, and estimation method. With the response simulated from normal distribution, both models estimated using either full or conditional likelihood give the expected results (Figures 1 and 2). However, with the response simulated from gamma distribution and the sampling fraction n/N small, the misspecified normal model estimated using full likelihood indicates a spurious association between the response and the covariate. This is because the missing data act as extra parameters, allowing the optimization procedure to obtain a better fit to the skewed data. The corresponding sampling distributions for β1 shown in Figure 3 have become bimodal. It should be noted that the initial values given to the optimization algorithm corresponded always to the “correct solution” of no covariate effect, thus enabling the algorithm to find the correct mode. However, bimodality in the sampling distribution does not necessarily indicate that the likelihood functions given a single realized dataset would be bimodal.Table 2
Maximum likelihood estimates of parametersβ of the outcome models Yi∣Zi,β~N(β0+β1Zi,β22) and Yi∣Zi,β~NPQM(β0+β1Zi,β2,β3,β4) and the population frequency γ of the binary covariate Zi (means over 1000 replications). Outcome variable was simulated from Yi~N(4,4) and Yi~Gamma(4,1).
Normal datanNormal modelNPQM modelβ̂0β̂1SE(β̂1)β̂2γ̂β̂0β̂1SE(β̂1)β̂2β̂3SE(β̂3)β̂4γ̂Full lik.10004.000.000.162.000.204.000.000.162.000.000.010.140.205004.000.010.221.990.204.000.010.221.990.000.010.140.202004.000.000.351.990.204.000.000.351.990.000.010.140.201004.000.010.501.990.204.000.010.501.990.000.020.140.20Cond. lik.5004.000.010.221.990.204.000.010.221.990.000.020.140.202004.010.000.351.990.204.010.000.351.990.000.030.130.201004.000.010.501.980.204.000.010.501.980.000.050.130.20Gamma datanNormal modelNPQM modelβ̂0β̂1SE(β̂1)β̂2γ̂β̂0β̂1SE(β̂1)β̂2β̂3SE(β̂3)β̂4γ̂Full lik.10004.000.010.161.990.203.990.000.132.010.150.010.140.205004.000.010.221.990.203.990.010.202.010.150.010.140.202003.970.150.381.980.203.980.070.352.000.150.020.140.201003.860.910.381.880.193.950.180.321.960.160.020.140.22Cond. lik.5004.000.000.221.990.204.000.000.192.010.150.020.140.202004.00-0.010.351.980.204.010.000.292.000.150.030.140.201004.00-0.010.501.970.204.010.010.412.000.150.050.130.20Figure 1
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; normal model fitted to normal data (Yi~N(4,4)).Figure 2
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; NPQM model fitted to normal data (Yi~N(4,4)).Figure 3
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; normal model fitted to gamma data (Yi~Gamma(4,1)).Using the NPQM model which allows a skewed residual distribution does not correct the situation either when combined with full likelihood estimation (Figure4). This is not unexpected since adding more parameters to an already overparameterized situation is not necessarily helpful in solving identification problems. In contrast, conditional likelihood continues to give reasonable results in all cases, with the skewness of the data correctly reflected by the parameter β3 of the NPQM model. Some degree of misspecification of the response model is required for the multimodality to appear since when we simulated the response variable from NPQM distribution, the full likelihood fit of the NPQM model indicated no problems, while the full likelihood fit of the normal model showed the same problems as with gamma data (Figures 5 and 6). We also repeated all the simulation alternatives with β1=1 (Yi~Zi+N(4,4) and Yi~Zi+Gamma(4,1)), with the conclusions essentially unchanged (Table 3 and Figures 7, 8, 9, and 10).Table 3
Maximum likelihood estimates of parametersβ of the outcome models Yi∣Zi,β~N(β0+β1Zi,β22) and Yi∣Zi,β~NPQM(β0+β1Zi,β2,β3,β4) and the population frequency γ of the binary covariate Zi (means over 1000 replications). Outcome variable was simulated from Yi∣Zi~Zi+N(4,4) and Yi∣Zi~Zi+Gamma(4,1).
Normal datanNormal modelNPQM modelβ̂0β̂1SE(β̂1)β̂2γ̂β̂0β̂1SE(β̂1)β̂2β̂3SE(β̂3)β̂4γ̂Full lik.10004.001.000.162.000.204.001.000.162.000.000.010.140.205004.001.010.221.990.204.001.010.221.990.000.020.140.202004.001.000.341.990.204.001.000.341.990.000.020.140.201004.001.000.471.990.204.001.010.481.990.000.020.140.20Cond. lik.5004.001.010.221.990.204.001.010.221.990.000.020.140.202004.011.000.351.990.204.011.000.351.990.000.030.130.201004.001.010.501.980.204.001.010.511.980.000.050.130.20Gamma datanNormal modelNPQM modelβ̂0β̂1SE(β̂1)β̂2γ̂β̂0β̂1SE(β̂1)β̂2β̂3SE(β̂3)β̂4γ̂Full lik.10004.001.010.161.990.203.991.000.132.010.150.010.140.205003.971.150.241.980.203.971.080.182.000.150.020.140.202003.832.010.361.860.193.921.280.231.970.160.020.140.211003.722.950.261.700.163.871.440.261.940.160.020.140.23Cond. lik.5004.001.000.221.990.204.001.000.192.010.150.020.140.202004.000.990.351.980.204.011.000.292.000.150.030.140.201004.000.990.501.970.204.011.010.412.000.150.050.130.20Figure 4
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; NPQM model fitted to gamma data (Yi~Gamma(4,1)).Figure 5
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; normal model fitted to NPQM data (Yi~NPQM(4,2,0.15,0.14)).Figure 6
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; NPQM model fitted to NPQM data (Yi~NPQM(4,2,0.15,0.14)).Figure 7
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; normal model fitted to normal data (Yi~Zi+N(4,4)).Figure 8
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; NPQM model fitted to normal data (Yi~Zi+N(4,4)).Figure 9
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; normal model fitted to gamma data (Yi~Zi+Gamma(4,1)).Figure 10
Sampling distributions for maximum likelihood estimates of regression coefficientβ1 from 1000 replications; NPQM model fitted to gamma data (Yi~Zi+Gamma(4,1)).


## 6.3. An Example with Real Data

The case-cohort set for all-cause mortality in the example cohort (N=5039) is of size n=1816, the union of a subcohort of size 1068 and 996 deaths due to any cause (since the subcohort has been selected independently of the case status, it includes a number of cases). The case-cohort design applied in this cohort was described in Section 4.3.1. Due to sample unavailability or unsuccessful genotyping, the genotype for the lactase persistence SNP rs4988235 was unavailable for 156 individuals, giving |𝒪∖ℳ|=1660. The model for the data described in Section 1 was specified as follows: Hardy-Weinberg model Zi∣γ~Binom(2,γ) was used for the number of lactase persistence alleles. The model for BMI, which is the response variable of interest, given age at baseline Xi≡bi and the genotype Zi, was taken to be Yi∣bi,Zi,β~NPQM(β0+β1bi+β21{Zi=0},β3,β4,β5). The model for the survival outcome was taken to be P(Ti∈dti,Ei=ei∣Yi,Zi,α)∝α[λi(ti)]1{ei=1}Si(ti), where Si(ti)≡exp{∫0tiλi(u)du}, and using again the Weibull form for the hazard function λi(u)≡(α4/α0)(u/α0)α4-1exp{α1bi+α2Yi+α31{Zi=0}}. In both regression models, the lactase persistence allele noncarriers are compared to the hetero- and homozygote carriers of the allele. Under the case-cohort design where the subcohort sampling probabilities depend on bi, the conditional likelihood expression takes the form   (6.2)P(TO∖M,EO∖M,YO∖M∣RC,MO,TC∖O∪M,EC∖O∪M,YC∖O∪M,XC,ZO∖M,θ)∝(α,β)∏i∈O∖MP(Ti,Ei∣Yi,Zi,α)P(Yi∣bi,Zi,β)∫yi[1-(1-π(bi))Si(ci)]P(dyi∣bi,Zi,β),
or   (6.3)P(MO,TO,EO,YO,ZO∖M∣RC,TC∖O,EC∖O,YC∖O,XC,θ)∝(α,β,γ)∏i∈O∖MP(Ti,Ei∣Yi,Zi,α)P(Yi∣bi,Zi,β)P(Zi=zi∣γ)∑zi∫yi[1-(1-π(bi))Si(ci)]P(dyi∣bi,zi,β)P(Zi=zi∣γ)×∏i∈M∑ziP(Ti,Ei∣Yi,zi,α)P(Yi∣bi,zi,β)P(Zi=zi∣γ)∑zi∫yi[1-(1-π(bi))Si(ci)]P(dyi∣bi,zi,β)P(Zi=zi∣γ),
depending on whether we condition upon the observed genotype data Z𝒪∖ℳ, or whether this is modeled as part of the likelihood (see Section 5.1.3). The difference between these two approaches is that the former requires more assumptions on the missingness mechanism. Therefore, we compare here the two approaches to see whether there is any observable difference between the results. In addition, expression (6.3) enables estimation of the population allele frequency γ. However, it should be noted that while the summation over the single SNP genotype variable Zi is not a computational problem in the present example, this is not necessarily the case generally, with multiple continuous covariates in the model. Therefore, the potential advantage of expression such as (6.2) is that it avoids specification of the covariate distribution and the integration over the missing covariate values.A remaining issue to be considered is the numerical evaluation of the integrals overyi in the correction terms of (6.2) and (6.3). It should be noted that these have to be evaluated at each parameter value tried at the hill-climbing numerical optimization. Although numerical (quadrature) integration would be feasible in one dimension, here we propose evaluation of such integrals using Monte Carlo integration by repeated sampling from the distribution P(Yi∈dyi∣bi,Zi,β). This has an added advantage that while the numerical evaluation of the density of the NPQM model is computationally expensive, random variates can be easily simulated from this distribution since it is defined through the quantile function. This proceeds by drawing random variates from [0,1]-uniform distribution and applying formula (6.2) of Karvanen [20]. With random variates yi(j)~NPQM(β0+β1bi+β21{Zi=0},β3,β4,β5), j=1,…,k, drawn for each i∈𝒪∖ℳ given the current parameter values, the integrals in (6.2) are then approximated by the means (1/k)∑j=1k[1-(1-π(bi))P(Ti>ci∣yi(j),Zi,α)].Maximum likelihood estimates obtained by numerical maximization of the expressions (6.2) and (6.3) with respect to α, β, and in the latter case also γ, are presented in Table 4 (α̂ and γ̂) and Table 5 (β̂). For comparison, we also maximised a full likelihood of the form (5.1) with respect to (α,β,γ) and a weighted pseudolikelihood of the form (5.2) with respect to β. In the latter case, we accounted for the missingness within the case-cohort set by fitting a logistic regression model logit{P(Mi=0∣Ei,bi,Yi,η)}=η0+η1Ei+η2bi+η3Yi in the set 𝒪 and calculated the adjusted weights as inverses of P(Mi=0∣Ei,bi,Yi,η̂)P(Ri=1∣Ei,bi). Due to the extra estimation step, we used bootstrap with 5000 replications to obtain standard errors for the weighted pseudolikelihood estimates. Both conditional likelihood expressions gave very similar results, suggesting that the missing data within the case-cohort set is not a major issue in the present case. The comparison between different number of Monte Carlo replicates in numerical evaluation of the conditional likelihood correction term suggests that the estimates do no longer appreciably change when k is increased from 1000. However, the change in the estimates when increasing k from 100 to 1000 is already well within the standard errors.Table 4
Maximum likelihood estimates (standard errors) of the parametersα of the all-cause mortality model λi(t)=(α4/α0)(t/α0)α4-1exp{α1bi+α2Yi+α31{Zi=0}}, and γ, the population frequency of the lactase persistence allele.
Likelihoodkα̂0logα̂0α̂1α̂2α̂3α̂4logα̂4γ̂Conditional10023.883.17 (0.21)0.12 (0.01)0.00 (0.01)-0.01 (0.12)1.340.29 (0.03)—(6.2)100023.523.16 (0.21)0.12 (0.01)0.00 (0.01)-0.01 (0.12)1.340.29 (0.03)—500023.453.15 (0.21)0.12 (0.01)0.00 (0.01)-0.01 (0.12)1.340.29 (0.03)—Conditional10020.433.02 (0.20)0.12 (0.01)-0.01 (0.01)-0.01 (0.11)1.330.29 (0.03)0.59 (0.01)(6.3)100020.343.01 (0.20)0.12 (0.01)-0.01 (0.01)-0.01 (0.11)1.330.29 (0.03)0.59 (0.01)500020.443.02 (0.20)0.12 (0.01)-0.01 (0.01)-0.01 (0.11)1.330.29 (0.03)0.59 (0.01)Full—18.402.91 (0.17)0.09 (0.01)-0.01 (0.01)-0.07 (0.11)1.330.28 (0.03)0.59 (0.01)Table 5
Maximum likelihood estimates (standard errors) of the parametersβ of the BMI model Yi∣bi,Zi,β~NPQM(β0+β1bi+β21{Zi=0},β3,β4,β5).
Likelihoodkβ̂0β̂1β̂2β̂3log(β̂3/π)β̂4β̂5logβ̂5Conditional10026.91 (0.14)-0.07 (0.02)-0.64 (0.27)4.270.88 (0.02)0.14 (0.03)0.37-1.00 (0.05)(6.2)100026.92 (0.14)-0.07 (0.02)-0.64 (0.27)4.270.88 (0.02)0.14 (0.03)0.37-1.00  (0.05)500026.92 (0.14)-0.07 (0.02)-0.64 (0.27)4.270.88 (0.02)0.14 (0.03)0.37-1.00 (0.05)Conditional10026.93 (0.13)-0.07 (0.02)-0.63 (0.26)4.250.87 (0.02)0.14 (0.03)0.37-1.00 (0.05)(6.3)100026.93 (0.13)-0.07 (0.02)-0.63 (0.26)4.250.87 (0.02)0.14 (0.03)0.37-1.00 (0.05)500026.93 (0.13)-0.07 (0.02)-0.63 (0.26)4.250.87 (0.02)0.14 (0.03)0.37-1.00 (0.05)Full—26.94 (0.07)-0.07 (0.01)-0.60 (0.25)4.060.83 (0.01)0.13 (0.01)0.35-1.06 (0.03)Weighted—26.94 (0.15)-0.06 (0.03)-0.68 (0.33)4.300.88 (0.04)0.17 (0.05)0.41-0.89 (0.11)Full likelihood and weighted pseudolikelihood estimates agreed well with the conditional likelihood ones, although the latter had higher standard errors, as was the case also in the simulations. As noted by Kettunen et al. [10], the absence of the lactase persistence allele shows association with lower body mass index. The residuals from the model for BMI are significantly skewed to the right, as indicated by the estimates of β4.


## 7. Discussion

Although conditional logistic likelihood is well known in the context of risk set sampling designs (e.g., [5, 18, 40]), its connection to the general concept of conditional likelihood has rarely been emphasized by authors. Furthermore, use of the conditional likelihood approach is not limited to the binary event outcome response situation. In this paper we have attempted to highlight these issues in the context of modeling a continuous secondary response variable under cohort sampling designs. Naturally, the conditional likelihood approach is also valid for the primary time-to-event outcome. Although here we considered only parametric specifications for the time-to-event outcome, semiparametric maximum likelihood techniques could also be utilized here. A potential disadvantage of the likelihood-based methods for secondary analysis is that they require specification of a model for the primary time-to-event outcome even though it is not of main interest in the secondary analysis setting, whereas specification of this model is avoided in the weighted pseudolikelihood approach.Compared to the full likelihood approach which would be applicable under any cohort sampling design, the advantage of the conditional likelihood is that modeling of the population distribution of the covariates collected in the second phase can be avoided, if this is not of primary interest. In addition, as we demonstrated in a simulation example, full likelihood expressions with most of the covariate data unobserved may no longer be well behaved, a problem which the conditional likelihood approach does not have. The disadvantage of the conditional likelihood approach is that the second-phase sampling mechanism, that is, the joint distribution of the inclusion indicators, needs to be specified. In the case-cohort/Bernoulli sampling situation, this is straightforward, but in mechanisms such as risk set sampling, where the sampling probabilities are specified only implicitly, resolving the joint sampling probability can be computationally nontrivial, and approximations may be needed in practice. This is avoided in the full likelihood approach, since it does not require specification of the sampling mechanism. The inverse-probability-weighted estimating function only requires specification of the first-order selection probabilities, but the possible dependencies induced by the sampling mechanism need to be addressed in the variance estimation step. The full and conditional likelihood approaches gave equivalent efficiencies in our simulated setting, although the result might be different if the first-phase data involves covariates which are highly predictive of the second-phase covariate of interest. In any case, both likelihood-based methods gave a clear improvement in efficiency in the secondary analysis setting compared to the inverse-probability-weighted method.



---

*Source: 931416-2012-04-24.xml*