Secondary Analysis under Cohort Sampling Designs Using Conditional Likelihood

1 Department of Epidemiology, Biostatistics and Occupational Health, McGill University, Montreal, QC, Canada H3A 1A2 2 Indic Society for Education and Development (INSEED), Nashik, Maharashtra 422 011, India 3 Department of Vaccines, National Institute for Health and Welfare, 00271 Helsinki, Finland 4 Department of Mathematics and Statistics, University of Tampere, 33014 Tampere, Finland 5 Department of Mathematics and Statistics, University of Helsinki, 00014 Helsinki, Finland


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 casecontrol designs.As extensions to the basic setting, we consider treatment of missing secondphase 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.

Notation, Assumptions, and Full Likelihood
To cover our motivating example and also the general case, we introduce first some notation.Let the set C ≡ {1, . . ., N} represent the individuals in the cohort.Primary outcome in the study is a time-to-event outcome characterized by random variables T i , E i , i ∈ C, where T i corresponds to event time and E i to the event type of individual i, with E i 0 indicating a censoring event and E i 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 Y i , here BMI, is observed on all study participants at the outset of the study.In addition, there may be other covariates X i , available on all i ∈ C, 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 Z i are collected only on the second-phase study group O ≡ {i : R i 1} ⊆ C, specified by the inclusion indicators R i ∈ {0, 1}, analogously to the survey response/nonresponse setting of Rubin 21 .We will henceforth use vector notations of the type Z O ≡ {Z i : i ∈ O} 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 C , T C , E C , X C , Y C , Z O .We are interested in the model P Y i | X i , Z i , β for the secondary outcome, more precisely, the parameter describing the association between Y i and Z i , 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 X i -covariates; then if all subsequent analyses are conditional on X i , 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 C .This is assumed to be unconfounded with What follows could be further generalized by assuming the sampling mechanism to be ignorable so that 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 T i , E i , Y i , X i , Z i 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 R i .Now, following, for example, Saarela et al. 8 , we can write a full observed data likelihood expression

2.1
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 C\O which may be large compared to O , as well as modeling of the population distribution of Z i .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.

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 probabilities 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 In a typical cohort sampling design, all cases would receive unit weights, while noncases selected to the set O 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.

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 indicators R C , which produced the second-phase study group O 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 into V or W at will, depending on the parameters of interest.For instance, if the parameters γ are not of interest, we may condition on Z O by placing it in V .If the parameters γ need to be estimated, Z O may be placed in W.
3 We must have 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, page 156 , 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.

Conditional Likelihood Expression: The General Case
Following the stated rules, and making the same general assumptions as in Section 2, we proceed by partitioning the observed data into relevant subsets and, using a shorthand notation of Q i ≡ T i , E i , Y i for all outcome variables, work with a conditional likelihood 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 R i .
The ratio of the numerator and the second term in the denominator can be further written as

4.2
Here, the product forms follow from the exchangeability assumption for the random vectors Q i , X i , Z i .Thus, we have obtained where the numerator factors into the parametric models 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 C , E C , Y C , X C .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 where θ is given by 4.2 , and we have

4.5
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 Z i , and it can be easily seen that if there is no association between T i , E i and Z i and between Z i and X i , we lose no information relevant to learning about the association between Y i and Z i , 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 C\O , 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 Z i .In Section 6.1, the efficiencies of expressions 2.1 and 4.3 are compared in a simulated setting.

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 C or both X C , Y C in the subcohort is approximately matched to that of the cases by first fitting a logistic regression model, say, logit{P E i 1 | X i , μ } μ X i , and then selecting the subcohort using Bernoulli sampling with the probabilities π i ≡ 1/ 1 exp{− μ X i } , independently of the case status.We then have , where in practice, we make the approximation μ 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∈C π i ≈ mπ i / NP E i 1 .More generally, the subcohort selection probability π i may be any known function 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 firstorder approximation when the sampling fractions are small.In the Bernoulli sampling case, the conditional likelihood correction term reduces into the product form Let E ≡ {i : E i 1} denote the set of individuals who died during the followup, and let d ≡ |E|.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 i ∈ E, while to increase the efficiency of the design, the subcohort is selected with probabilities which depend on the age b i at the start of the followup included in the X i covariates through a logistic model as discussed above, giving the inclusion probability for non-cases as Under Bernoulli sampling, the realized sample size n ≡ |O| is random, with the expected sample size in the present example given by In the case of a mortality outcome and type I censoring at predetermined times c i , with the observed time given by T i ≡ min T i , c i , where T i is the latent time of death, 4.6 further simplifies into see also 19, pages 12-13

4.7
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 casecohort design E i 1 is not all-cause mortality, but, say, mortality due to cardiovascular diseases.Deaths due to other causes, denoted as E i 2, then appear as censoring.As before, let T i denote the latent time of death, with E i indicating the type of death or end of followup period.Similarly as above, we get

4.8
Here, the probabilities , are given by the cumulative incidence functions and P T i ≥ c i | X i , Y i , Z i , α 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 O, 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.

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, and m controls per case j are selected using without replacement sampling from the risk set R j ≡ {i : where A i T j is the at-risk indicator for subject i at event time T j .This is carried out independently for all cases j ∈ E. Let the sampled set of time-matched controls for case j be denoted as S j ⊆ R 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 S E ≡ {S j : j ∈ E} denote the collection of the sampled risk sets and S ≡ j∈E S j the pooled set of sampled controls.Noting that knowing all of T C , E C specifies the risk sets R E ≡ {R j : j ∈ E} as well as the order of the events, the risk set sampling mechanism is specified by the joint probability distribution 4.9 10

Journal of Probability and Statistics
If the sampling has been carried out within strata specified by the covariates Y C and/or X C through some function g X i , Y i , the risk sets would be defined within strata as R j {i : A i T j 1, g X i , Y i g X j , Y j } \ {j}, and the above reasoning would apply similarly with the redefined risk sets.The conditional likelihood correction term now becomes 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 O. 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 ∈ C \ E given by Thus, it might be tempting to replace 4.9 However, the properties of such an approximation will require further research.

Missing Second-Phase Covariate Data
Further issue to be considered in practical applications is possible missing covariate data within the set O. 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 O ≡ {M i : i ∈ O} be a set of indicator variables indicating whether the measurement of Z i turned out to be unsuccessful after selection of subject i into the set O, and let M ≡ {i : M i 1} ⊆ O. Now Z i is observed only on the set O \ M. We consider the implications of this separately for all of the three approaches discussed above.

Full Likelihood
If the missingness can be assumed missing at random, that is, the missingness indicators can again be included in the proportionality constant, and the full likelihood expression becomes

Inverse Probability Weighting
If the missingness mechanism would be known and missing at random , a weighted pseudolikelihood expression for the set O\M could be obtained by multiplying the original selection probabilities by the probabilities of observing the value of Z i after selection In practice, however, the missingness mechanism would have to be modeled to obtain estimates for these probabilities.

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 C and O 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 secondphase 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 Sections 1 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 T i denotes the observed time of the first incident disease event, death, or censoring, with E i 0 indicating censoring and E i ∈ {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.

Full Likelihood
The full likelihood expression 2.1 is now conditioned on the selection rule T C ≥ b C ≡ T i ≥ b i for all i ∈ C, that is, event-free survival until the age at the cohort baseline.The likelihood expression becomes

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 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.

Conditional Likelihood
As in Section 4.2, we get where the correction term becomes

5.7
In the case-cohort design discussed in Section 4.3.1, with type I censoring, this further simplifies into 5.8

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 Section 1 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 C , E C , X C , Y C and predetermined associations with the primary and secondary outcomes, an additional "missing" covariate Z C is simulated from the posterior predictive distribution Case-cohort sampling was then simulated in the complete datasets T C , E C , X C , Y C , Z C 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 Y i available and eligible for the case-cohort study.The primary time-toevent outcome T i , E i was taken to be all-cause mortality, with d 996.Given the primary and secondary outcome data and age at baseline X i ≡ b i observed for the cohort, and the current parameter values, additional binary covariate values were drawn from the conditional distributions where the model for the survival outcome, 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 model 3 was used for the secondary outcome, while the population distribution of the additional covariate was specified as 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 Y i and b i were centered but otherwise untransformed.1000 values for each Z i , i ∈ C, 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 Z i , 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 1761.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 Even though this is a different model than the one specified for the simulation step, we did not simulate the timeto-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 y i 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.

Multimodality under Full Likelihood When the Sampling Fraction Is Small
Let now the observed data be only R C , Y C , Z O , and let the sampling mechanism be simple random sampling without replacement; P R C 1/ N n , 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 C \ O, to misspecification of the response model.The full likelihood expression now becomes Due to the simple random sampling, the corresponding conditional likelihoods simplify into 1000, we simulated covariate values from Z i ∼ Bernoulli 0.2 and two different sets of response values both independently of Z i from Y i ∼ N 4, 4 and Y i ∼ 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 Z i | γ ∼ Bernoulli γ , with the models for the response specified as 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 Conditional likelihood, n = 100 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.
Table 2 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.Conditional likelihood, n = 100 Using the NPQM model which allows a skewed residual distribution does not correct the situation either when combined with full likelihood estimation Figure 4 .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 Y i ∼ Z i N 4, 4 and Y i ∼ Z i Gamma 4, 1 , with the conclusions essentially unchanged Table 3 and Figures 7, 8, 9, and 10 .

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.        Conditional likelihood, n = 100 sample unavailability or unsuccessful genotyping, the genotype for the lactase persistence SNP rs4988235 was unavailable for 156 individuals, giving |O \ M| 1660.The model for the data described in Section 1 was specified as follows: Hardy-Weinberg model Z i | γ ∼ 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 X i ≡ b i and the genotype Z i , was taken to be

Table 3: Maximum likelihood estimates of parameters β of the outcome models
The model for the survival outcome was taken to be 0 λ i u du}, and using again the Weibull form for the hazard function 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 b i , the conditional likelihood expression takes the form Conditional likelihood, n = 100 or depending on whether we condition upon the observed genotype data Z O\M , 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 Z i is not a computational problem in the present example, this is not necessarily the case generally, with multiple continuous covariates in the model.Therefore, Conditional likelihood, n = 100 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 over y i 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 Y i ∈ dy i | b i , Z i , β .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 y j i 1, . . ., k, drawn for each i ∈ O \ M given the current parameter values, the integrals in 6.2 are then approximated by the means 1/k 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 Conditional likelihood, n = 100 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.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 .

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  Conditional likelihood, n = 100 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 timeto-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 Conditional likelihood, n = 100 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.Conditional likelihood, n = 100 be the correct data generating mechanism.This becomes Conditional likelihood, n = 100 is correct, that is, corresponds to the data generating mechanism.Thus, the expression 3.1 gives a valid estimating equation for parameters β.However, since the variance of the pseudoscore function does not equal the Fisher information, a robust/sandwich type variance estimator must be used instead of the inverse of the observed information.As an example, below we derive this in the special case of a Bernoulli sampling mechanism , where this is straightforward.For consideration of asymptotic variances of inverse-probability-weighted estimators under stratified without replacement sampling or risk set sampling, we refer to Breslow and Wellner 41 and Cai and Zheng 42 .
Denoting q β ≡ i∈O log P Y i | X i , Z i , β /P R i | T i , E i , Y i , X i and β ≡ arg max β q β , the latter is given as a solution to the estimating equation q β 0, or approximately, by the first-order Taylor expansion at the true value β 0 , as a solution to q β 0 q β 0 β − β 0 0, which in turn, by substituting expected to observed pseudoinformation, gives the approximate relationship β − β 0 ≈ E −q β 0 −1 q β 0 .Taking covariance of both sides then gives cov β ≈ E −q β 0 −1 cov q β 0 E −q β 0 −1 , where the pseudo-Fisher information is given by Since the pseudoscore function has zero mean and now R i ⊥ R j , i / j, its covariance is given by cov q β 0 E ⎡ ⎣ i∈C j∈C In practice, these quantities are evaluated at the maximum likelihood point, replacing the expected quantities by their observed counterparts as Alternatively, bootstrap variance estimation might be utilized, but it should be noted that standard bootstrap would be valid only under Bernoulli type sampling designs, which can be interpreted as sampling from infinite population, whereas dependencies induced by sampling without replacement would require application of finite population bootstrap methods e.g., 43, 44, page 37 .

B.1
The last form is the retrospective form suggested by Lin and Zeng 7, page 257 for secondary analysis under case-control designs.Thus, retrospective likelihood is a special case of conditional likelihood.The advantage of the above expression is that the terms specifying the sampling mechanism cancel out.The drawbacks are that the baseline risk level part of α is not identifiable from this expression see, e.g., 30, page 1076 and that the parameters γ have to be estimated or otherwise made to disappear Lin and Zeng 7 applied profile likelihood .However, nothing prevents us from choosing a more useful partitioning of the observed data.

B.3
The baseline risk level may be identifiable from this last expression, depending on the sampling mechanism used cf.18 .The misconception that case-control design necessitates the use of retrospective likelihood has been recently criticized by Langholz 45 and results from equating likelihood function to a data generating mechanism.However, the order of the data collection does not need to determine how the likelihood is factorized, as long as the sampling mechanism is properly taken into account.

C. Mean and Variance of the Conditional Likelihood Score Function
In the following, we suppress the covariates X C from the notation since these are always conditioned upon and do not affect the derivations.Also, for notational simplicity, the parameter θ is taken to be a scalar, and the notation is compressed so that, for example, P dz O | θ ≡ P Z O ∈ dz O | θ ≡ P Z i ∈ dz i for all i ∈ O | θ .The probability measures indexed with respect to θ are assumed to possess the regularity properties allowing interchanging the order of integration with respect to the random variables and differentiation with respect to θ.
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 β 3 L
-skewness , and λ 4 β 4 L-kurtosis .This distribution is suitable for regression modeling as the first L-moment is the mean of Full likelihood,

Figure 1 :
Figure 1: Sampling distributions for maximum likelihood estimates of regression coefficient β 1 from 1000 replications; normal model fitted to normal data Y i ∼ N 4, 4 .

Figure 2 :
Figure 2: Sampling distributions for maximum likelihood estimates of regression coefficient β 1 from 1000 replications; NPQM model fitted to normal data Y i ∼ N 4, 4 .

4 γ
4 and the population frequency γ of the binary covariate Z i means over 1000 replications .Outcome variable was simulated from Y Full lik.

4 and 4 γ
the population frequency γ of the binary covariate Z i means over 1000 replications .Outcome variable was simulated from Y Full lik.

Figure 3 :
Figure 3: Sampling distributions for maximum likelihood estimates of regression coefficient β 1 from 1000 replications; normal model fitted to gamma data Y i ∼ Gamma 4, 1 .

Figure 4 :
Figure 4: Sampling distributions for maximum likelihood estimates of regression coefficient β 1 from 1000 replications; NPQM model fitted to gamma data Y i ∼ Gamma 4, 1 .

Figure 7 :
Figure 7: Sampling distributions for maximum likelihood estimates of regression coefficient β 1 from 1000 replications; normal model fitted to normal data Y i ∼ Z i N 4, 4 .

Figure 8 :
Figure 8: Sampling distributions for maximum likelihood estimates of regression coefficient β 1 from 1000 replications; NPQM model fitted to normal data Y i ∼ Z i N 4, 4 .

Figure 9 :
Figure 9: Sampling distributions for maximum likelihood estimates of regression coefficient β 1 from 1000 replications; normal model fitted to gamma data Y i ∼ Z i Gamma 4, 1 .

Full
last form being the expectation of the complete data score function for parameters β, which becomes zero through the inner expectation if the parametric model P Y i | X i , Z i , β

Figure 10 :
Figure 10: Sampling distributions for maximum likelihood estimates of regression coefficient β 1 from 1000 replications; NPQM model fitted to gamma data Y i ∼ Z i Gamma 4, 1 .
Consider a special case where the observed data are only R C , E C , Y O , Z O , where E i , i ∈ C, are binary indicators for the case status.This corresponds to the case-control situation considered by Jiang et al. 6 and Lin and Zeng 7 .The set C now represents a study base from which the cases and controls have been selected.Now choosing a partitioning W Y O , Z O and V R C , E C , and supposing only the case status information has been used in the selection of cases and controls, we obtain a conditional likelihood C\O , X C and model the population distribution of Z i 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.
θ is obtained from 4.5 by replacing the set O with O \ M. 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 M are either missing completely at random or that M O depends only on X C , 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 C , the missingness mechanism would have to be modeled.However, in such a case, it may be easier to work with the partitioning WM O , T O , E O , Y O , Z O\M and V R C , T C\O , E C\O , Y

Table 1 :
Maximum likelihood estimates of the parameters α of the all-cause mortality model

Table 2 :
Maximum likelihood estimates of parameters β of the outcome models Y

Table 4 :
Maximum likelihood estimates standard errors of the parameters α of the all-cause mortality model λ

Table 5 :
Maximum likelihood estimates standard errors of the parameters β of the BMI model Y For instance, with W Y O and V R C , E C , Z O , we obtainP Y O | R C , E C , Z O , θ P R C | E C P E C , Y O , Z O | θ P R C | E C P E C , Z O | θ i∈O P E i | Y i , Z i , α P Y i | Z i , β y i P E i | y i , Z i , α P Y i ∈ dy i | Z i , β O , Y O and V R C , E C\O , Z O , analogously to Section 4.2, P E O , Y O | R C , E C\O , Z O , θ C , Y O , Z O | θ P R C | E C\O , θ P E C\O , Z O | θ i∈O P E i | Y i , Z i , α P Y i | Z i ,β y i :i∈O e i :i∈O P R C | E C i∈O P E i e i | y i , Z i , α P Y i ∈ dy i | Z i , β .
θ ∝ P E The conditional likelihood function of interest is nowP dq O | r C , q C\O , z O , θ P r C | q C ,z O , θ P dq C , dz O | θ P r C | q C\O , z O , θ P dq C\O , dz O | θ P dq O , dq C\O | z O , θ P dz O | θ P r C | q C\O , z O , θ P dq C\O | z O , θ P dz O | θ P dq O| q C\O , z O , θ P r C | q C\O , z O , θ , C.1and the corresponding score function becomes∂ log P dq O | r C , q C\O , z O , θ ∂θ ∂ log P dq O | q C\O , z O , θ ∂θ − ∂ log P r C | q C\O , z O , θ ∂θ ∂P dq O | q C\O , z O , θ /∂θ P dq O | q C\O , z O , θ − ∂P r C | q C\O , z O , θ /∂θ P r C | q C\O , z O , θ .The expectation of the score function is considered with respect to the full model P r C | q C , z C , θ P dq C , dz C | θ .The expectation of the first term in C.2 becomesE ∂P dq O | q C\O , z O , θ /∂θ P dq O | q C\O , z O , θ r C q C z C ∂P dq O | q C\O , z O , θ /∂θ P dq O | q C\O , z O , θ P r C | q C , z C , θ P dq C , dz C | θ r C q C z O ∂P dq O | q C\O , z O , θ /∂θ P dq O | q C\O , z O , θ × z C\O P r C | q C , z C\O , z O , θ P dq C , dz C\O , dz O | θ r C q C z O ∂P dq O | q C\O , z O , θ /∂θ P dq O | q C\O , z O , θ × P r C | q C , z O , θ P dq O | q C\O , z O , θ P dq C\O , dz O | θ .In the above expressions, it should be noted that the index set O is fixed by r C .Due to the assumption on unconfounded sampling mechanism, the term P r C | q C , z O , θ does not depend on θ, and we getE ∂P dq O | q C\O , z O , θ /∂θ P dq O | q C\O , z O , θ r C q C z O ∂ P r C | q C , z O , θ P dq O | q C\O , z O , θ ∂θ P dq C\O , dz O | θ r C q C\O z O ∂ q O P r C | q O , q C\O , z O , θ P dq O | q C\O , z O , θ ∂θ P dq C\O , dz O | θ r C q C\O z O ∂P r C | q C\O , z O , θ ∂θ P dq C\O , dz O | θ .C\O , z O , θ P r C | q C\O , z O , θ P dq C\O , dz O | θ r C q C\O z O ∂P r C | q C\O , z O , θ ∂θ P dq C\O , dz O | θ q O P dq O | r C , q C\O , z O , θ r C q C\O z O ∂P r C | q C\O , z O , θ ∂θ P dq C\O , dz O | θ .C\O , z O , θ , C.6and the conditional likelihood score function has zero expectation.The Fisher information is given byE − ∂ 2 log P dq O | r C , q C\O , z O , θ ∂θ 2 E ∂ 2 log P r C | q C\O , z O , θ ∂θ 2 − E ∂ 2 log P dq O | q C\O , z O , θ ∂θ 2 E ∂ 2 P r C | q C\O , z O , θ /∂θ 2 P r C | q C\O , z O , θ − E ∂P r C | q C\O , z O , θ /∂θ P r C | q C\O , z O , θ 2 E ∂P dq O | q C\O , z O , θ /∂θ P dq O | q C\O , z O , θ P dq O | q C\O , z O , θ /∂θ 2 P dq O | q C\O , z O , θ .By writing open the expectations as above, it is easy to see thatE ∂ 2 P dq O | q C\O , z O , θ /∂θ 2 P dq O | q C\O , z O , θ E ∂ 2 P r C | q C\O , z O , θ /∂θ 2 P r C | q C\O , z O , θ , ∂ 2 log P dq O | r C , q C\O , z O , θ ∂θ 2 E ∂P dq O | q C\O , z O , θ /∂θ P dq O | q C\O , z O , θ C | q C\O , z O , θ /∂θ P r C | q C\O , z O , θ 2 E ∂P dq O | q C\O , z O , θ /∂θ P dq O | q C\O , z O , θ C | q C\O , z O , θ /∂θ P r C | q C\O , z O , θ C | q C\O ,z O , θ /∂θ P r C | q C\O , z O , θ ∂P dq O | q C\O , z O , θ /∂θ P dq O | q C\O , z O , θ E ∂P dq O | q C\O , z O , θ /∂θ P dq O | q C\O , z O , θ − ∂P r C | q C\O , z O , θ /∂θ P r C | q C\O , z O , θ θ ∝ C | q C\O , z O , θ /∂θ P r C | q C\O , z O , θ r C q C z C ∂P r C | q C\O , z O , θ /∂θ P r C | q C\O , z O , θ P r C | q C , z C , θ P dq C , dz C | θ r C q C z O ∂P r C | q C\O , z O , θ /∂θ P r C | q C\O , z O , θ × P dq O | r C , q O | q C\O , z O , θ /∂θ P dq O | q C\O , z O , θ E ∂P r C | q C\O , z O , θ /∂θ P r C | q E ∂P r C | q C\O , z O , θ /∂θ P r C | q C\O , z O , θ ∂P dq O | q C\O , z O , θ /∂θ P dq O | q C\O , z O , θ r C q C\O z O ∂P r C | q C\O , z O , θ /∂θ P r C | q C\O , z O , θ ∂P r C | q C\O , z O , θ ∂θ P dq C\O , dz O | θ E ∂P r C | q C\O , z O , θ /∂θ P r C | q C\O , z O , θ 2 .C.8Thus, the Fisher information becomes E − that is, equal to the variance of the conditional likelihood score function.