Bayesian subcohort selection for longitudinal covariate measurements in follow‐up studies

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Bayesian subcohort selection for longitudinal covariate measurements in follow‐up studies © 2022 The Authors. Statistica Neerlandica published by John Wiley & Sons Ltd on behalf of Netherlands Society for Statistics and Operations Research. Published version Reinikainen, Jaakko; Karvanen, Juha Reinikainen, J., & Karvanen, J. (2022). Bayesian subcohort selection for longitudinal covariate measurements in follow‐up studies. Statistica Neerlandica, 76(4), 372-390. https://doi.org/10.1111/stan.12264 2022 Received: 25 April 2019 Revised: 20 December 2021 Accepted: 17 January 2022 DOI: 10.1111/stan.12264 ORIGINAL ARTICLE Bayesian subcohort selection for longitudinal covariate measurements in follow-up studies Jaakko Reinikainen1,2 Juha Karvanen2 1Finnish Institute for Health and Welfare, Helsinki, Finland 2Department of Mathematics and Statistics, University of Jyvaskyla, Jyvaskyla, Finland Correspondence Jaakko Reinikainen, Finnish Institution forHealthandWelfare,Helsinki,Finland. Email: jaakk[email protected] Funding information Emil Aaltosen Säätiö Abstract We propose an approach for the planning of longitudinal covariate measurements in follow-up studies where covariates are time-varying. We assume that the entire cohort cannot be selected for longitudinal measurements due to financial limitations, and study how a subset of the cohort should be selected optimally, in order to obtain precise estimates of covariate effects in a survival model. In our approach, the study will be designed sequentially utilizing the data collected in previousmeasurementsoftheindividualsaspriorinforma- tion.WeproposeusingaBayesianoptimalitycriterionin the subcohort selections, which is compared with simplerandomsamplingusingsimulatedandrealfollow-up data. Our work improves the computational approach compared to the previous research on the topic so that designswithseveralcovariatesandmeasurementpoints can be implemented. As an example we derive the optimaldesignforstudyingtheeffectofbodymassindexand smokingonall-causemortalityinaFinnishlongitudinal study. Our results support the conclusion that the precisionof theestimates can be clearly improvedby optimal design. KEYWORDS Bayesian optimal design, data collection, follow-up study, longitudinal measurements, study design This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2022 The Authors. Statistica Neerlandica published by John Wiley & Sons Ltd on behalf of Netherlands Society for Statistics and Operations Research. 372 wileyonlinelibrary.com/journal/stan Statistica Neerlandica. 2022;76:372–390. REINIKAINEN and KARVANEN 373 1INTRODUCTION Longitudinalcovariatemeasurementsareoftencarried outin follow-up studieswhen thecovariates are time varying. These measurements give useful information about the trajectories of the covariates. Frequent remeasurements provide more information than infrequent ones, but in practice, limited resources may restrict the number of measurements and researchers have to consider how to design the study cost efficiently. We study optimal design in a scenario where we cannot afford to remeasure the entire cohort but can only select a subset of the cohort, called a subcohort. The goal is to estimate the effects of the covariates on survival as precisely as possible. The study is designed sequentially, which here means that the subcohorts are selected just before the measurement times and all information collected prior to a new measurement is utilized. This kind of design procedure is realizable if researchers can clearly define beforehand the purpose of data collection, that is, the parameters of interest to be estimated from the data. In addition, the outcome data must already be available during the follow-up. The proposed method requires especially that up-to-date survival information can be obtained when needed. This is possible, for example, in Finland, where data on mortality and hospitalizations are available from administrative registries. Use of a Bayesian version of Ds-optimality, an optimality criterion based on Fisher information, is proposed for the selection of the subcohort. In addition to the Ds-criterion, there are also other optimality criteria, which were originally developed for design of experiments (Atkinson, Donev, & Tobias, 2007; Pukelsheim, 1993), but can also be applied in observational studies. For example, Karvanen, Kulathinal, and Gasbarra (2009) considered optimal subset selection for genotyping in a follow-up study, Buzoianu and Kadane (2009) investigated selection of patients for a diagnostic test and Mehtälä, Auranen, and Kulathinal (2015) studied optimal time spacings for observations of a multistate Markov process. Reinikainen, Karvanen, and Tolonen (2016) studied the same problem with a frequentist approach that was restricted to one time-varying covariate and one remeasurement after the baseline. The limitation of this approach is that there is no straightforward way to generalize it to more realistic scenarios where several covariates and measurement points are allowed. The Bayesian approach presented in this paper offers flexibility that makes it possible to handle these challenges. In a multivariate setting with multiple longitudinal measurements, several computational challenges are faced. Because of the very large discrete design space, it is practically impossible to evaluate all possible designs. Thus, some heuristic method is needed to find an approximate solution.BayesiananalysiswithMarkovchainMonteCarlo(MCMC)offersflexibilityinthemodeling,butrequirescomputationalresources.Selectingonlysubcohortsforremeasurementcreates a large amount of missing data, which further complicates the data analysis. Thedatamissingbydesignhavepreviously(Reinikainenetal.,2016)beenhandledwithmultiple imputation and a likelihood-based approach with numerical integration. These approaches do not look promising for the generalized problem because multiple imputation of covariates conditioned on survival data would be complicated and numerical integration would become infeasible as increasing the number of covariates increases the dimension of the integral. Here, we handle missing data by using Bayesian data augmentation, which is expected to be a more flexible method when the number of measurement points and covariates is increased. Parallel computing is utilized in a simulation study. 374 REINIKAINEN and KARVANEN The optimal subcohort selection is studied using simulated and real data. The real data consist of the Finnish cohorts from the Seven Countries Study (Keys, 1970), an international epidemiologic study characterized by a long follow-up time and several longitudinal covariate measurements. We use body mass index and smoking as time-varying covariates and all-cause mortality as the outcome. With these data, the proposed selection procedure is compared with simple random sampling (SRS) of individuals to be remeasured and with a case where the entire cohort is selected for remeasurement. 2SURVIVAL MODEL Inthissection,weintroducethenotationforourstudy designandsurvivalmodel.Thesearelater used to present our optimal subcohort selection procedure. Many of the following assumptions are made to be suitable for the real data example of Section 6, but the general idea is applicable to other designs and models as well. Let us consider a follow-up study in which survival time is the response variable and M longitudinal measurements are carried out for the time-varying covariates after the baseline measurement. We denote the covariate values by xmjh for the measurement m=0,…,M,the individual j=1,…,Nand the covariate h=1,…,H. The corresponding random variables are denoted by Xmjh and for all the covariates and individuals shortly by Xm.The measurements for the covariates are carried out at time points 𝜏0,…,𝜏Min calendar time for the individuals who arealiveandhavebeenselectedtobe measured.Theindividualsinthecohortmaybeofdifferent ages at the time of the baseline measurement and the start of the follow-up 𝜏0. The follow-up has a predetermined length ending at the time 𝜏M+1. The exact survival times are available from a registry. The observed survival information at the time of mth remeasurement for the individual jis denoted by ymj =(tmj,𝛿mj),where tmj is the continuously measured survival time with age as the time scale and 𝛿mj is the status indicator tellingthe vitalstatus attime 𝜏m(𝛿mj =1 for an eventand 𝛿mj =0 for censoring).Thus, if𝛿mj =1, tmj is the age at the event and if 𝛿mj =0, tmj is the age at 𝜏m. Although the survival time is observed continuously, the piece-wise modeling approach uses separate time and status variables for each measurement time interval. The individual jhas time-to-event variables t1j,t2j,…for each part of the follow-up where they are still alive. We denote the random variables related to survival information by Ymj and for all the individuals by Ym.We use notation xm=(xm1,…,xmH)Tand 𝜷=(𝛽1,…,𝛽H)Tand assume that the covariate effects fulfill the Markov property 𝜆(tm+1|x0,…,xm,𝛿m=0)=𝜆(tm+1|xm,𝛿m=0). We continue by assuming that the covariates are related to the hazard of the event through the proportional hazards model 𝜆(tm+1|xm,𝛿m=0)=𝜆0(tm+1|𝛿m=0)exp(𝜷Txm).(1) Conditioning on 𝛿m=0 means that only those individuals contribute here who have not died before the mth remeasurement. As the covariate information is updated during the follow-up, also the hazard of an individual is assumed to change correspondingly. The survival times are REINIKAINEN and KARVANEN 375 assumed to follow the Weibull distribution, when the baseline hazard function has the form 𝜆0(tm+1|𝛿m=0)=a b(tm+1 b)a−1, whereaistheshapeparameterandbisthescaleparameter.Thenmodel(1)becomesaparametric form of the time-dependent Cox model (Therneau & Grambsch, 2000). Other distributions than the Weibull could also be used for the survival times. Now, we can write the survival function and the density function as S(tm+1|xm,𝛿m=0)=S0(tm+1|𝛿m=0)exp(𝜷Txm)and p(tm+1|xm,𝛿m=0)=𝜆(tm+1|xm,𝛿m=0)S(tm+1|xm,𝛿m=0), where S0(tm+1|𝛿m=0)is the baseline survival function. The modeling is carried out piecewisely in time, because the covariate information changes at the measurement points. For this reason we have to deal with left-truncated and possibly right-censored Weibull distributions in the time intervals (𝜏0,𝜏1],(𝜏1,𝜏2],….Survival times are left-truncated at the lower limit of a time interval byscalingthelikelihoodwithS(tmj|xmj,𝛿m=0)becauseanobservedsurvivaltimetm+1,jcannotbe smallerthantmj.Thisisa sequential conditioningwherethelikelihoodisscaledwith the survival probability at age tmj.In addition, left-truncation is necessary because the age of an individual is used as the time scale and we do not assume that the follow-up would begin at time zero, which would be the time of birth. Thus, we have to take into account that persons who have died before thestartof the follow-updidnothave anopportunitytobeincluded inthecohort.The likelihood contribution for the individual jfor the parameters 𝛽1,…,𝛽H,aand bis Lj(𝛽1,…,𝛽H,a,b)= m′ j ∏ m=0(p(tm+1,j|xmj,𝛿m=0) S(tmj|xmj,𝛿m=0))𝛿m+1,j(S(tm+1,j|xmj,𝛿m=0) S(tmj|xmj,𝛿m=0))1−𝛿m+1,j ,(2) where m′ j=max{0,…,M∶𝛿mj =0}. 3OPTIMAL SUBCOHORT SELECTION If the entire cohort cannot be remeasured because of financial limitations, we have to select a subcohort, which we can afford to measure. The optimal selection aims to make the estimates of the parameters of interest as precise as possible subject to financial constraints. In this paper, we focus on the estimation of regression parameters 𝜷. Assume that baseline covariate measurements (and possibly some longitudinal measurements) have been carried out and that continuous survival information can be obtained during the study for all individuals. When we want to carry out the next longitudinal measurement for a subcohort, we proceed by taking the following general steps: 1. Just before the new longitudinal measurement, use the data already collected to obtain prior information about the parameters of interest 2. Define the optimality criterion as an expectation over the prior distribution 3. Maximize the optimality criterion and select an optimal subcohort for the remeasurement 4. Remeasure the covariates for the selected subcohort 376 REINIKAINEN and KARVANEN The rest of this section and Section 4.1 provides a description of one possible way to perform steps 2. and 3. Section 4.2 focuses on step 1. 3.1 Bayesian optimal design in a general form An optimal design problem can be seen as a problem of maximizing the expected utility U(𝝃) for a design 𝝃from a design space Ξ(Chaloner & Verdinelli, 1995). See the seminal paper by Lindley(1956)onBayesianoptimaldesignandWoods,Overstall,Adamou,andWaite(2017)fora recent review on the topic. In our problem, 𝝃is an indicator matrix with individuals on rows and measurement times on columns, where an element (j,m)is 1 if individual jhas been selected for measurementmand0otherwise.Theconstraintoflimitedresourcesmeansherethatthecolumn sums are fixed in 𝝃.The column sums need not be the same, because we may have different amount of resources for different re-examinations. Data w=(x,y),wherexcorresponds to covariate data and yto survival outcome, come from a sample space .The outcome data yand the baseline covariate measurements are assumed to be available on all individuals, whereas longitudinal covariate data is collected according to the design 𝝃.The data are assumed to follow a model p(w|𝜽),where parameters 𝜽belong to the parameter space Θ. The fully Bayesian solution for the optimal design problem would involve integrating a measureofobservedutilityoverdatawandtheposteriordistributionofparameters𝜽.Insteadofthis, we use a common approach (Atkinson et al., 2007; Chaloner & Verdinelli, 1995) where the integration is done over the prior density distribution of 𝜽and the utility is defined as a function of the expected information. With this notation, the expected utility can be written as U(𝝃)=∫Θg[Ew|𝜽,𝝃{Iw(𝜽)}]p(𝜽)d𝜽,(3) where gis a function such as determinant (D-optimality), Ew|𝜽,𝝃{Iw(𝜽)} is an expectation of the Fisher infromation of the parameters 𝜽over data wand p(𝜽)is the prior density distribution of 𝜽. Here,weusetheDs-criterion(Atkinsonetal.,2007),whichisaspecialcaseofthewidelyused D-optimality criterion. The D-optimal design maximizes the determinant of the Fisher information matrix or equivalently minimizes the determinant of the covariance matrix. Ds-optimality considersonly a subset of sparameters.If theparametervector𝜽=(𝜃1,…,𝜃s,…,𝜃p)Tincludes first the sparameters of interest and then p−snuisance parameters, Ds-optimal design minimizes the determinant of the s×supper left submatrix of I(𝜽)−1.In our case, sis the number of 𝛽-parameters in the survival model, and thus we will call the criterion D𝛽-optimality. 3.2 Selection criterion for the Weibull proportional hazards model Let us consider the Fisher information matrix of the model introduced in Section 2including parameters 𝜽∗=(𝛽1,…,𝛽H,a,b): IX,Y(𝜽∗)=−E(𝜕2logp(X0,…,XM,Y1,…,YM+1) 𝜕𝜽∗2),(4) REINIKAINEN and KARVANEN 377 where Ym=(Tm,𝛿m). The parameter vector 𝜽∗includes the parameters of the survival model but not the parameters of covariate processes. We assume Markov properties p(Ym+1|X0,…,Xm,Y1,…,Ym)=p(Ym+1|Xm,Ym)and p(Xm|X0,…,Xm−1,Y1,…,Ym)=p(Xm|Xm−1), and decompose the logarithmic joint distribution as logp(X0,…,XM,Y1,…,YM+1)=logp(X0)+logp(Y1|X0)+…+logp(XM|XM−1) +logp(YM+1|XM,YM),which allows us to decompose IX,Y(𝜽∗)similarly. The Markov assumption is made here to avoid unnecessarily complicating the description of the method. Other dependencystructurescouldbeassumedaswell,dependingonthedataforwhichthemethodisapplied. When we are selecting individuals for the mth remeasurement, we have measured covariates X0,…,Xm−1, which may include missing values and survival information is observed up to the timeofthe mth remeasurement, that is,Ymisknown. Therefore,theselectionforthe mth remeasurementisbasedontheexpectationsofXmandYm+1utilizingthepreviouslyobserveddata.Only those individuals who have not yet died or been censored can be considered as candidates for the remeasurement. That is, optimization is carried out by selecting individuals who are still in the risksetforfurtherexamination.Individualswhosecovariateshavenotbeenmeasuredpreviously are also candidates for the remeasurement. Due to the Markov assumptions, it is sufficient that the selection is based only on the expectations with respect to the next unobserved part and not of all the forthcoming parts of the follow-up. Now, using the above-mentioned decomposition in matrix (4), the information matrix used in the selection for the mth remeasurement can be written as Im X,Y(𝜽∗)=−E[𝜕2 𝜕𝜽∗2{logp(X0)+logp(Y1|X0)+… +logp(Xm−1|Xm−2)+logp(Ym|Xm−1,Ym−1)} +E{(𝜕2 𝜕𝜽∗2logp(Xm|Xm−1)+ 𝜕2 𝜕𝜽∗2logp(Ym+1|Xm,Ym))|||X0,…,Xm−1,Ym}]. In the presentation of this matrix, the last term includes unobserved values, which is why an additionalconditionalexpectationistakenonit.Above,thetermswhichdonotincludevariableY vanish, because they do not include the survival model parameters 𝜽∗(parameters of the Weibull proportional hazards model). This leads to Im X,Y(𝜽∗)=E[−𝜕2 𝜕𝜽∗2{logp(Y1|X0)+···+logp(Ym|Xm−1,Ym−1)}] +E[E{−𝜕2 𝜕𝜽∗2logp(Ym+1|Xm,Ym)|||X0,…,Xm−1,Ym}] =IY1|X0(𝜽∗)+···+IYm|Xm−1,Ym−1(𝜽∗)+E{IYm+1|Xm,Ym(𝜽∗)},(5) where the outer expectation of the last term is with respect to unobserved data Ym+1|Xm,Ym. As values Y1,…,Ymand X0,…,Xm−1are already observed, the first mterms in (5)are replacedbytheobservedinformationJ(𝜽∗).Theinformationmatrixisthenamixtureofobserved and expected information and its element in row iand column kis Ψm X,Y(𝜽∗)i,k=JY1|X0(𝜽∗)i,k+···+JYm|Xm−1,Ym−1(𝜽∗)i,k+E{IYm+1|Xm,Ym(𝜽∗)}i,k =−N ∑ j=1[𝜕2 𝜕𝜃∗ i𝜕𝜃∗ k logp(y1j|x0j)] 378 REINIKAINEN and KARVANEN −···− Nm−1 ∑ j=1[𝜕2 𝜕𝜃∗ i𝜕𝜃∗ k logp(ymj|xm−1,j,ym−1,j)] − nm ∑ j=1[E{𝜕2 𝜕𝜃∗ i𝜕𝜃∗ k logp(Ym+1,j|Xmj,ymj)|||x0,j,…,xm−1,j,ymj}],(6) where Nm−1is the number of individuals who have not had an event or been censored before the measurement m−1andnmis the number of individuals to be selected for the mth measurement. Note that when only a subcohort has been selected for the measurements at a time point 𝜏m′,m′∈{1,…,m−1},thencovariatesaremissingfortheindividualsnotselected.Wedescribe in Section 4.2, how these missing data are handled. The subcohort selection is carried out just before the new measurement, so Ym+1and Xmare not observed for anyone. The expectation can be calculated by Monte Carlo integration. The calculation of Ψm X,Y(𝜽∗)requires the second-order partial derivatives of logp(y1|x0),…,logp(ym+1|xm). The value of the D𝛽-criterion is obtained by taking the determinant of the H×Hupper left submatrix of Ψm X,Y(𝜽∗)−1,where 𝜽∗=(𝛽1,…,𝛽H,a,b).We denote this value of the criterion by Dm 𝛽(𝝃m,𝜽∗),where 𝝃mis an indicator matrix describing which individuals have been measured at the time points 𝜏0,…,𝜏mand 𝜽∗emphasizes that the criterion depends on the parameters. 3.3 Bayesian selection Now, we combine the Bayesian optimal design theory introduced in Section 3.1 and the criterion derivedin Section 3.2. Here,theoptimalselectionispresentedin a generalformandthepractical solution for the search problem is given in Section 4.1. Fisher information matrices of nonlinear models usually depend on model parameters (Chaloner & Verdinelli, 1995), which is also the case in our application. Therefore, some prior information about the parameters is needed in order to use the D𝛽-criterion. In the Bayesian approach the information obtained from the data already collected during the follow-up and/or from previous studies can be used to provide prior distributions of the parameters when applying the optimality criterion. In other words, we use informativepriors,whichareactuallyposteriorsfromthedataalreadycollected,inthesubcohort selections. In the selection for the mth remeasurement, the prior probability density distribution pm(𝜽∗)is, in fact, the posterior p(𝜽∗|x0,…,xm−1,y1,…,ym). To minimize the D𝛽-criterion, we specify the last column of 𝝃mso that the expected utility U(𝝃m)=−∫𝜽∗Dm 𝛽(𝝃m,𝜽∗)pm(𝜽∗)d𝜽∗,(7) willbemaximized. Above, Dm 𝛽(𝝃m,𝜽∗)isthevalueoftheD𝛽-criterionforthe mth remeasurement depending on the design 𝝃mand the parameters 𝜽∗. This is a specific form of the Equation (3). 4COMPUTATIONAL IMPLEMENTATION 4.1 Search for optimal design Theintegralin(7)canbeapproximatedbysamplingparametervaluesfromthemultivariateprior probabilitydensitydistributionpm(𝜽∗),generatingdata(Xm,Ym+1)giventheparametersandthen REINIKAINEN and KARVANEN 379 replacingtheintegralwithamean 1 q∑q l=1Dm 𝛽(𝝃,𝜽∗ l),whereqisthenumberofrealizationssampled from pm(𝜽∗)(Atkinson, Demetrio, & Zocchi, 1995). The priors become more informative during the follow-up as the amount of collected data increases. In addition to the model parameters and the new data, the missing covariate values are also treatedasunknownparameters.Thepredictivedistributionsofthemissingvaluesareusedasthe prior distributions in the selections. This means that if the previous measurements include missing data, the criterion (7) averages also over these informative prior distributions of the missing values. We draw qrealizationsfromthe priorsofthe missing valuesand use themsimilarlytothe realizations from pm(𝜽∗). In practice, the number of different subcohorts that could be selected for the new measurementiseasilysolargethatitiscomputationallyimpossibletogothrougheachofthem.Therefore some heuristic method is needed. We use a so-called greedy method (Wright & Bailer, 2006), also known as sequential search (Dykstra, 1971), to find an approximately optimal subcohort. This method selects nindividuals sequentially one by one: when k−1 individuals (0<k<n)have been selected for the subcohort, the kth selection is made so that the expected utility is maximizedonthecondition thatinformationfrompreviouslyselected k−1 individualsisincludedin the calculation of the criterion. The procedure goes on similarly by selecting the next individual to be included in the subcohort so that the expected utility is maximized taking into account the information obtained from the previously selected individuals. From the beginning of the selection procedure, the information matrix (6) includes all the informationalreadycollectedduringthefollow-up.Alltheindividualswhohavenothadanevent areconsideredascandidatesforthenewmeasurement.Theprocedurecontinuesbytestingwhich candidate should be included in the expectation part on the last row of (6) in order to obtain the minimum value of the Bayesian optimality criterion Dm 𝛽(𝝃m,𝜽∗). If there are two or more individuals who would minimize the criterion, the selection between them can be done randomly. The expected information of the selected individual is then included permanently into (6)and the procedure continues until the subcohort has reached the predetermined size. The covariate measurements are carried out after the whole subcohort has been selected. 4.2 Parameter estimation The estimation of the parameters 𝜽∗of the survival model (2) is needed before each subcohort selectionandfinallywhenthefollow-upstudyhasended.Itshouldbeemphasizedthatthemodel used in the subcohort selections need not be the same as the model used in the final analysis of the data, although the subcohort will be optimal with respect to the selection model. Once the data have been collected, the validity of the model assumptions can be reassessed and the model can be changed if needed. If the variables are the same in both models, the selected subcohort is likely tobe betterthan a simple randomsample, althoughthe finalanalysis modelcould bemore complex than the selection model. This is exemplified in Section 6. Inthesubcohortselections,theestimatedposteriordistributionsoftheparametersareusedas informative prior density distributions. The parameters are estimated using a Bayesian approach with MCMC sampling. We use noninformative priors for the parameters: 𝜷∼N(0,104⋅I),r∼ Gamma(1,0.0001)and 𝛼∼N(0,104),where rand 𝛼are reparameterized Weibull parameters, so that r=aand 𝛼=−alogb,0is a zero vector and Iis an H×Hidentity matrix. The subcohort selection may lead to a large amount of data “missing by design,” which we handle by Bayesian data augmentation (Tanner & Wong, 1987). The data missing by design are 386 REINIKAINEN and KARVANEN 0 100 300 500 20 25 30 35 Smokers Selection order 2nd meas. of BMI 0 100 300 500 20 25 30 35 Non−smokers Selection order 2nd meas. of BMI 0 100 300 500 20 25 30 35 2nd meas. of BMI missing Selection order Mean of the imputations of the 2nd meas. of BMI FIGURE 3 Selection order of individuals for the third measurement using D𝛽-optimality and the East–West data. These panels represent the right panel of Figure 2decomposed into smokers, nonsmokers and those who have missing value of body mass index (BMI) in the second measurement. Individuals with missing BMI are indicated by triangle symbols regardless of the missingness of the smoking status. Usually, the missingness of BMI means also the missingness of the smoking status in the data TABLE 3 Results for comparisons of using different models in the subcohort selections and in the final analysis using the East–West data. 500 individuals were selected for the second and third measurements. For simple random sampling (SRS), 𝛽1,𝛽2and SE are means of the posterior means and standard errors estimated from the Markov chain Monte Carlo chains from 1,000 analyses BMI (linear) BMI (quadratic) Analysis model Selection 𝜷1(SE) 𝜷2(SE) Quadratic Full cohort −0.059 (0.0082) 0.0073 (0.0011) SRS −0.059 (0.0096) 0.0079 (0.0013) D𝛽(quadr. model) −0.059 (0.0085) 0.0077 (0.0011) D𝛽(lin. model) −0.059 (0.0087) 0.0073 (0.0011) Linear Full cohort −0.038 (0.0077) SRS −0.033 (0.0084) D𝛽(quadr. model) −0.038 (0.0084) D𝛽(lin. model) −0.037 (0.0084) Abbreviation: BMI, body mass index. revealsthat therearelesssmokersthannonsmokersin the selectedsubcohort. Infact,therewere 104smokers,203nonsmokersand459individualswithmissingsmokingstatusamongthecandi- dates for the third measurement, of which all 104 smokers, 138 nonsmokers and 358 individuals with missing smoking status were selected for the third measurement. Table 2shows that there is clear benefit of using the D𝛽-design instead of the SRS. All the standard errors in the D𝛽-design are smaller than in the SRS-design for each subcohort size used. The D𝛽-selection leads usually to estimates closer to those obtained using the full cohort, than the SRS. Surprisingly, both selection methods seem to lead to greater estimates of the effect of smoking than the full cohort. REINIKAINEN and KARVANEN 387 The standard errors of the quadratic term of BMI and smoking, obtained using SRS, are smallerwhenn=300thanwhenn=400,whichisanunexpectedresult.Threehundredindividualsareonly24%oftheindividualsaliveatthetimeofthesecondmeasurementand 39% of those alive at the time of the third measurement, which may be too small proportions in a real study to obtain reliable estimates. Estimation may become sensitive to model misspecification when the proportion of missing data becomes large (Saarela, Kulathinal, & Karvanen, 2012). In practice, an analyst would not necessarily like to use the same model in optimal selection and in the final analysis. Table 3shows the results from an example of using different models in theselectionsand final analysis.Here,weusethesamedataasinthepreviousexample,but(centered) BMI as the only covariate. When the analysis model is quadratic, the quadratic selection model leads to slightly better precision than the linear selection model but both selection models still outperform SRS. When the analysis model is linear, the results do not deteriorate even if the quadratic selection model is used. 7DISCUSSION The cost-efficiency of a follow-up study can be improved by careful planning of longitudinal measurements. The present paper considered the scenario where we can afford measuring the time-varyingcovariatesonlyforasubsetofthecohort.WeproposedusingaBayesianapproachin optimalsubcohortselectionwithaFisherinformationbasedD𝛽-optimality.Ourworkovercomes thelimitationsofthepreviouspaper(Reinikainenetal.,2016),whereasimplescenariowithonly one covariate and one re-measurement was considered. The estimates and their precisioncorresponding to the D𝛽-selection and SRS were compared. The use of the D𝛽-optimality led to more precise estimates and the precision was seen to remain satisfactory compared with the full cohort design. The results indicated that in order to obtain estimates as precise as possible for regression parameters of the survival model, old individuals with extreme covariate values should be preferred, which is consistent with our previous results (Reinikainen et al., 2016). A similar result was obtained when we used a covariate with quadratic effect (BMI) in the real follow-up data example.Anothercovariateusedinthisexamplewassmokingstatusasabinaryvariable.Optimal selection seemed to balance the expected number of smokers and nonsmokers measured. The general idea of measuring only a subcohort is applied in many epidemiological study designs, like case-control and case-cohort designs and their variants (Keogh & White, 2013; Kulathinal et al., 2007; Sun, Joffe, Chen, & Brunelli, 2010). However, in our application the setup is different. We use an approach with an explicit utility function describing the goal of the study. The presented approach borrows elements from optimal design of experiments and applies them to the design of an observational study. We recommend the use of a Bayesian approach in this kind of sequential study design problem. Some prior knowledge is always required in nonlinear design problems, because optimal designs depend on model parameters and so a Bayesian approach with informative priors is a natural way to incorporate this knowledge into design optimization. Bayesian data augmentation appears here to be a flexible method for the handling of missing data when we increase the number of measurement points and covariates. The proposed approach can be adopted to other scenarios where repeated measurements are conducted.However,itisdifficulttoprovidea generalpurposesoftwareimplementationbecause the details of Bayesian models differ case by case and the forms of survival model and optimality 388 REINIKAINEN and KARVANEN criterion should be selected to fit the available data and the aims of the study. We recommend that researchers test the performance of their implementations in simulations similar to those presented in Table 1before applying them in real study design. Inexperiments,the optimal designforalinear model mightconsistofonlytwodesign points, which would give no power for detecting nonlinear effects. In principle, the same applies also in our setup, but in practice, the problem is not realized in observational studies with continuous covariates and moderate sample size. The reason for this is that only a few individuals with the optimal covariate values are available in the cohort and after they are selected, the selection procedure must rely on individuals with a wider variety of covariate values. Although the covariate processes were assumed to be piece-wise constant in the subcohort selections, a joint model of survival and longitudinal data (Rizopoulos, 2012) could be considered in the final analysis for more realistic treatment of time-varying covariates. In the general case with multiple longitudinal covariates, joint modeling is, however, computationally very demanding and would be further complicated in our approach with a large amount of missing data. We considered the design optimization with respect to only one model at a time. However, if two or more models or utility functions would be of interest, compound design criteria could be applied (Atkinson et al., 2007). Then, the optimality criterion should include the parameters of all models of interest. If we later want to use the collected data to some purpose not addressed in subcohort selection, the optimality does not hold anymore, and in an extreme case the selected subcohort could perform even worse than SRS. This situation is unlikely to occur, if a new outcome variable has the same covariates as the one used in optimal selection, or if new covariates are correlated with those used in optimization. The situation is similar to case-control studies where the controls for a specific outcome can be used for other outcomes (Saarela et al., 2012; Saarela, Kulathinal, Arjas, & Läärä, 2008) although they are not optimal. The selection procedure presented in this paper requires that up-to-date survival information is already available during the study. The information of the measured covariates is also needed during the study if it is used in the optimal selection. The proposed approach is also applicable to some retrospective designs. For instance, consider a study where blood samples or other biological specimen are collected and stored for all individuals and years later some biomarkers are measured from the stored sample. Selecting only a subcohort for these measurements may be a reasonable option if the extraction of the biomarkers is expensive. Then an approach similar to one presented in this paper could be used to optimally select the subcohort. Theproposedselectionmethod is nonrandomin the sensethat individuals areselected deterministically according to the selection criterion, but we do not see this as a disadvantage when theparametersof a survivalmodelareof interest.Distributionsof the covariates orabsoluterisks of the outcome in the population can also be assessed using the estimated analysis model. If it is important to assess the distributions of the covariates without relying on model assumptions, the subcohort should be selected as a random sample (Kulathinal et al., 2007). In this setting a randomized version of the sequential design construction could be considered (Atkinson & Biswas,2014). Our method mayintroducesome selectionbias, if the effect of a covariatechanges with age and this has not been taken into account in the selection. The procedure could also be developed to include the costs of data collection in the optimization (Fedorov & Leonov, 2013). ACKNOWLEDGEMENTS The research of the first author was supported by the Emil Aaltonen Foundation. Hanna Tolonen from Finnish Institute for Health and Welfare is acknowledged for providing the data of the REINIKAINEN and KARVANEN 389 East–West study. The authors thank CSC – IT Center for Science Ltd. for providing computing resources and Marie Reilly, Kari Auranen and Juho Kopra for helpful comments. DATA AVAILABILITY STATEMENT ThedatathatsupportthefindingsofthisstudyareavailableonrequestfromtheFinnishInstitute for Health and Welfare. The data are not publicly available due to privacy or ethical restrictions. ORCID Jaakko Reinikainen https://orcid.org/0000-0002-7825-5359 REFERENCES Atkinson, A., Demetrio, C., & Zocchi, S. (1995). Optimum dose levels when males and females differ in response. Journal of the Royal Statistical Society Series C (Applied Statistics),44, 213–226. Atkinson, A. C., & Biswas, A. (2014). Randomised response-adaptive designs in clinical trials. Boca Raton, FL: CRC Press. Atkinson,A.C.,Donev,A.N.,&Tobias,R.D.(2007).Optimumexperimentaldesigns,withSAS.Oxford,UK:Oxford University Press. Buzoianu, M., & Kadane, J. B. (2009). Optimal Bayesian design for patient selection in a clinical study. Biometrics, 65, 953–961. Chaloner, K., & Verdinelli, I. (1995). Bayesian experimental design: A review. Statistical Science,10, 273–304. Corrada, M. M.,Kawas, C. H., Mozaffar, F., & Paganini-Hill, A. (2006). Association of body mass index and weight change with all-cause mortality in the elderly. American Journal of Epidemiology,163, 938–949. Dykstra, O. (1971). The augmentation of experimental data to maximize ∣X′X∣.Technometrics,13, 682–688. Elfving, G. (1952). Optimum allocation in linear regression theory. The Annals of Mathematical Statistics,23, 255–262. Fedorov, V. V., & Leonov, S. L. (2013). Optimal design for nonlinear response models. Boca Raton, FL: CRC Press. Karvanen,J.,Kulathinal,S.,&Gasbarra,D.(2009).Optimaldesignstoselectindividualsforgenotypingconditional onobservedbinaryor survivaloutcomesand non-geneticcovariates.Computational Statistics &DataAnalysis, 53, 1782–1793. Keogh, R. H., & White, I. R. (2013). Using full-cohort data in nested case–control and case–cohort studies by multiple imputation. Statistics in Medicine,32, 4021–4043. Keys, A. (1970). Coronary heart disease in seven countries. Circulation,41, 186–195. Kulathinal,S.,&Arjas,E.(2006).Bayesianinferencefromcase-cohortdatawithmultipleend-points.Scandinavian Journal of Statistics,33, 25–36. Kulathinal, S., Karvanen, J., Saarela, O., Kuulasmaa, K., & MORGAM Project. (2007). Case-cohort design in practice – Experiences from the MORGAM project. Epidemiological Perspectives & Innovations,4, 15. Lindley, D. V. (1956). On a measure of the information provided by an experiment. The Annals of Mathematical Statistics,27, 986–1005. Lunn, D., Jackson, C., Best, N., Thomas, A., & Spiegelhalter, D. (2012). The BUGS book: A practical introduction to Bayesian analysis. Boca Raton, FL: CRC Press. Lunn, D., Spiegelhalter, D., Thomas, A., & Best, N. (2009). The BUGS project: Evolution, critique and future directions. Statistics in Medicine,28, 3049–3067. Mehtälä, J., Auranen, K., & Kulathinal, S. (2015). Optimal observation times for multistate Markov models–applicationstopneumococcalcolonizationstudies.JournaloftheRoyalStatisticalSociety:SeriesC(Applied Statistics),64, 451–468. Pukelsheim, F. (1993). Optimal design of experiments. New York, NY: Wiley. R Core Team. (2014). R: A language and environment for statistical computing. Vienna, Austria: Retrieved from: R Foundation for Statistical Computing. http://www.R-project.org/ Reinikainen, J., Karvanen, J., & Tolonen, H. (2016). Optimal selection of individuals for repeated covariate measurements in follow-up studies. Statistical Methods in Medical Research,25, 2420–2433. 390 REINIKAINEN and KARVANEN Reinikainen, J., Laatikainen, T., Karvanen, J., & Tolonen, H. (2015). Lifetime cumulative risk factors explain cardiovascular disease mortality in a 50-year follow-up study in Finland. International Journal of Epidemiology, 44, 108–116. Rizopoulos, D. (2012). Jointmodels forlongitudinal andtime-to-event data:With applicationsinR. Boca Raton, FL: CRC Press. Saarela, O., Kulathinal, S., Arjas, E., & Läärä, E. (2008). Nested case–control data utilized for multiple outcomes: A likelihood approach and alternatives. Statistics in Medicine,27, 5991–6008. Saarela, O., Kulathinal, S., & Karvanen, J. (2012). Secondary analysis under cohort sampling designs using conditional likelihood. Journal of Probability and Statistics,2012, 1–37. Sun,W.,Joffe,M.M.,Chen,J.,&Brunelli,S.M.(2010).Designandanalysisofmultipleeventscase–controlstudies. Biometrics,66, 1220–1229. Tanner, M. A., & Wong, W. H. (1987). The calculation of posterior distributions by data augmentation. Journal of the American Statistical Association,82, 528–540. Therneau, T. M., & Grambsch, P. M. (2000). Modeling survival data: Extending the Cox model. New York, NY: Springer. Woods, D. C., Overstall, A. M., Adamou, M., & Waite, T. W. (2017). Bayesian design of experiments for generalizedlinearmodelsand dimensionalanalysiswith industrialandscientific application.QualityEngineering,29, 91–103. Wright, S. E., & Bailer, A. J. (2006). Optimal experimental design for a nonlinear response in enviromental toxicology. Biometrics,62, 886–892. Zhao, W., Katzmarzyk, P. T., Horswell, R., Wang, Y., Li, W., Johnson, J., …Hu, G. (2014). Body mass index and the risk of all-cause mortality among patients with Type 2 diabetes. Circulation,130, 2143–2151. How to cite this article: Reinikainen, J., & Karvanen, J. (2022). Bayesian subcohort selection for longitudinal covariate measurements in follow-up studies. Statistica Neerlandica,76(4), 372–390. https://doi.org/10.1111/stan.12264