Second-Order Estimating Equations for Clustered Current Status Data from Family Studies Using Response-Dependent Sampling.

Studies about the genetic basis for disease are routinely conducted through family studies under response-dependent sampling in which affected individuals called probands are sampled from a disease registry, and their respective family members (non-probands) are recruited for study. The extent to which the dependence in some feature of the disease process (e.g., presence, age of onset, severity) varies according to the kinship of individuals reflects the evidence of a genetic cause for disease. When the probands are selected from a disease registry, it is common for them to provide quite detailed information regarding their disease history, but non-probands often simply provide their disease status at the time of contact. We develop conditional second-order estimating equations for studying the nature and extent of within-family dependence which recognizes the biased sampling scheme employed in family studies and the current status data provided by the non-probands. Simulation studies are carried out to evaluate the finite sample performance of different estimating functions and to quantify the empirical relative efficiency of the various methods. Sensitivity to model misspecification is also explored. An application to a motivating psoriatic arthritis family study is given for illustration.

Introduction

Introduction

The heritable nature of disease can be inferred from the structure of the within-family dependence in disease manifestation [19]. For rare diseases, population-based cohort studies are inefficient and impractical, so response-dependent biased family designs are routinely employed to obtain enriched samples with higher representation of diseased individuals and more variation in genetic markers than would be seen in the unselected population. Much work has been carried out for the analysis of such data when the disease status is modeled as a binary trait. Conditional likelihood based on generalized linear mixed models can be used for dealing with the dependence of binary phenotypes within families [3, 4], or estimating equations can be formed by specifying marginal mean and dependence structures for the analysis of binary phenotypes from case–control family studies [32].

The age of onset for many chronic diseases is highly variable, however, and simply using the binary trait of the disease status does not account for the variable times individuals have been at risk for disease in family studies. MacLean et al. [21] and Shih and Chatterjee [27] pointed out that, if information on the age of onset and the effect of censoring are not addressed, the estimators of the covariate effects on the disease process may be less efficient and the degree of familial aggregation may be underestimated. Models which consider the disease onset time distribution and measure dependence in terms of these times offer a preferable framework for analysis.

When interest lies in examining genetic association or gene-environment interaction, case–control or case-only family study are commonly used. Li et al. [20], Hsu et al. [15], and Shih and Chatterjee [27] proposed likelihood methods based on disease onset time for case–control family study and Chatterjee et al. [6] proposed methods to estimate the relative risk, cumulative risk, and residual familial aggregation for case–control family data and modified method for case-only family data. In their methods, modeling and estimation of the residual familial aggregation is key to adjustment for ascertainment bias, but this is done using an exchangeable dependence structure in which the association is the same for different pairs of relatives. Gorfine et al. [14] use the frailty models to account for heterogeneity in familial risk, but pointed out that frailty-based methods may be affected by the uncertainty on the frailty parameter estimate.

In this article, we consider a simple family study, where an affected individual called a proband is selected from a registry of patients. Consenting family members (non-probands) of each proband are then recruited and examined to collect information on their disease status [5]. Probands are given a special designation because their disease status led to the selection of their family. Under such sampling scheme, one obtains a right-truncated onset time for the proband and current status (type I interval-censored) data for non-probands [28]. While work has been done on the analysis of multivariate current status data [8, 16], little has been done to our knowledge in the context of biased sampling schemes.

Insights into the genetic basis of the disease can be gained by comparing the strength of the association in disease status between pairs of family members with different kinships [7, 18]. More elaborate dependence modeling also plays a central role when studying the “parent-of-descent” hypothesis, where the primary goal is to estimate and compare the strength of father–child and mother–child associations in phenotype to elucidate the role of the sex chromosomes in disease transmission [30]. With this in mind, we consider copula models [22] as a basis for modeling the joint risk of disease among family members. The dependence parameters can be interpreted as reflecting “residual familial aggregation” that is not explained by covariates in the marginal models. Copula models have several advantages over frailty models. First, the marginal models still retain simple interpretation when using copula models, which is not the case under the frailty model. Second, copula models yield dependence measures which are functionally independent of the parameters in the marginal onset time distribution, so the marginal distribution can be specified in any desirable way. Third, the dependence measure is directly specified under the copula model which has clear meaning and it also provides a natural basis for regression of genetic effects, but the frailty models do not provide simple measures of within-family dependence and it is difficult to interpret the meaning of the dependence.

Analyses must address the biased sampling scheme employed in these studies. Likelihood contributions from each family which are proportional to joint probability functions for the phenotypes of non-probands conditional on the disease status of the proband will admit valid inference [29] under correct model specification, but enumeration of all possible sample outcomes can be computationally demanding with large families. We develop a class of conditional second-order estimating equations in the spirit of Prentice [25]. We use the term conditional to reflect the fact that moments in the second-order estimating equation are all conditional on the disease onset time of the proband. A supplementary estimating equation is incorporated to extract limited information about the marginal onset time distribution from the proband.

The University of Toronto Psoriatic Arthritis Family Study

The incidence of psoriatic arthritis (PsA) is reported to be between 0.3 and 1.0% [9] and hereditary factors are thought to be important, as some studies have suggested that close blood relatives of individuals affected by psoriatic arthritis are at higher risk of developing the disease compared to the general population. Characterizing the within-family association nature and identifying important genetic risk factors are important to understand the disease etiology. Particular interest lies in assessing whether there is a higher rate of paternal, rather than maternal, transmission of the disease, which is also called “parent- of-origin” effect [2]. A family study of psoriatic arthritis is conducted in the Centre for Prognosis Studies in the Rheumatic Disease at the University of Toronto. Probands were selected from the members of the University of Toronto Psoriatic Arthritis Registry, and their family members were recruited into the family study with their consent. A total of 169 two-generation families ranging in size from 2 to 7 individuals were recruited; 54 families were comprised of only one non-proband and 115 have more than one non-proband. The disease onset times were only available for probands, but for other family members only the disease status is available when they are examined, yielding current status data. In total 538 individuals are in the family study and only 194 (169 probands and 25 non-probands) were diagnosed with PsA. Except for the demographic data, information of some HLA markers is also available for individuals in the PsA family study. We focus on identifying the significant HLA markers for the psoriatic arthritis and characterizing the within-family association structure, also testing whether there is “parent-of-origin” effect for the psoriatic arthritis.

The remainder of this paper is organized as follows. In Sect. 2 we define notation and formulate the conditional second-order estimating equation for family data under response-dependent sampling, which are a combination of right-truncated onset time from probands and current status data from non-probands. We consider an illustrative example in which the dependence structure is governed by a Gaussian copula and work with this model in subsequent calculations and simulations where we examine specific estimating equations involving different derivative matrices and working independence assumptions. In Sect. 3, we explore the asymptotic relative efficiencies and finite sample properties of estimators from several variants of the estimating equations introduced in Sect. 2; these results also permit sample size calculations for planning studies aiming to detect effects of genetic markers. The impact of misspecification of the dependence structure on properties of estimators and power of genetic tests is investigated in Sect. 4. An application to the motivating psoriatic arthritis family study is given in Sect. 5 in which we assess the genetic basis of the disease. Concluding remarks are given in Sect. 6.

Conditional Estimating Equations Under Biased Sampling

Notation, Sampling, and Observation Scheme for Family Studies

We consider the setting in which a registry of M individuals is created by selecting a random sample from a population, screening each individual for disease, and recruiting those found to have the condition of interest [10]. If denotes the age of individual 0 in family i at the time of sampling and screening, and denotes their age of disease onset, then this individual is recruited to the registry if ; we assume that is verifiable by a review of medical records for individuals recruited to the registry. When a family study is carried out, we assume that probands are selected from the disease registry by simple random sampling and without loss of generality we label the families of selected probands .

We let and denote, respectively, the event time and a covariate vector of individual j in family i, where are the labels for the non-probands. Then if and , we write the joint cumulative distribution function (j.c.d.f) for family i as . We assume , where , and write . The marginal hazard function for the disease onset time of individual j, , in family i iswhere we write under a proportional hazards formulation. This gives , where , is a vector, is a vector of regression coefficients, and . We let parameterize the within-family dependence and .

Classification of non-probands with respect to their disease status is made at the time of recruitment and clinical examination, yielding current status data. Let denote the age of non-proband j in family i at the time of assessment and let ; we let , and . If , , and , the family data therefore consist of subject to .

Second-Order Estimating Functions

The association parameter is of central importance here so we next formulate conditional second-order generalized estimating equations in the spirit of Prentice [25] and Zhao and Prentice [33].

Let be an vector of pairwise products of the elements in , where ; we let denote the element of corresponding to the pair (j, k) in family i. To account for response-biased sampling, we define conditional moments and let and be the contributions from the non-probands and let for the proband where we suppress the dependence on and . The conditional second-order estimating equations (CGEE2) denoted by have the form:withwhere , , , and , , , and ; note that unlike standard GEE2, since is functionally dependent on . The covariance matrices can be parameterized by the marginal and association parameters where the latter may be specified in terms of Kendall’s ; an example is given in Sect. 2.3. Consistent estimation of is possible based on the first term in (1), but the second term , where and , improves efficiency by exploiting the data on the onset time from the proband.

Subject to correct specification of the conditional moments, (1) is an unbiased estimating function, so the estimator solving is consistent with an asymptotic normal distributionwhereNatural empirical estimates of these matrices areandwhich yield .

Simplified forms of can be obtained by setting (denoted by ) or by letting both and (denoted by ). It is also common to simplify and adopt a form in which and while retaining the full structure of ; we refer to this as a working partial independence (WPI) matrix. Combining these simplifications, we consider four different estimating functions based on (1)

Full and Full covariance matrix denoted G–W,

Full and WPI denoted G-WPI,

G and WPI denoted G-WPI, and

G and WPI denoted -WPI.

An Illustrative Dependence Structure Based on a Gaussian Copula

The specific form of the moments for can be motivated by a copula model. Consider an vector of uniform [0, 1] variables , in which , . The j.c.d.f. for , denoted by , is a copula function in dimensions indexed by an parameter which characterizes the dependence [17, 22]. The Gaussian copula is a member of elliptical family of the form:where is the inverse cumulative distribution function of a standard normal random variable (r.v.) and is the j.c.d.f. of an multivariate normal r.v. with mean zero and correlation matrix ; is indexed by and we denote the off-diagonal entries by , . Specification of the Gaussian copula for induces a joint distribution for given bywhere , is a realization, and , . Copula functions such as this are attractive for dependence modeling since pairwise associations are parameterized to be functionally independent of the marginal parameters and different pairwise associations are permitted. The Kendall’s characterizing the association between and given , for example, is given by , . Regression modeling of the within-family dependence can be achieved by specifying a second-order model of the form , where is a 1–1 differentiable link function mapping Kendall’s onto the real line, is an covariate vector characterizing individuals j and k in family i and their relationship, and is the corresponding vector of coefficients. This second-order regression model can be helpful when investigating the effect of risk factors on the pairwise association as could represent family-level or individual-level features, or information on the kinship of individuals j and k in family i; inference on their effects can be easily carried out based on . For example, in the PsA family study with two generations, when the “parent-of-origin” hypothesis is of interest, we can formulate the second-order model asthen comparing and (or testing ) can inform us whether there is “parent-of-origin” effect in the onset time of PsA. More elaborate models which incorporate genetic covariates into the dependence model can also be specified.

Returning to the estimating function in (1), based on the Gaussian copula we have and , where , , and . The function is the j.c.d.f of a bivariate normal r.v. with mean zero and covariance matrix , whereThe entries of can also be derived based on the Gaussian copula where, for example, for withNote that this is a j.c.d.f for a multivariate normal r.v. with mean and covariance matrix denoted by , where andThese conditional moments are easily derived under a Gaussian copula.

Relative Efficiency Under Particular Estimating Equations

A Study of Asymptotic Relative Efficiency

Here we examine the asymptotic relative efficiency of four different conditional estimating equations as a function of the strength of the within-family association through the functions:where denotes an asymptotic variance and its subscript indexes the adopted conditional estimating equations proposed in Sect. 2.2 . All three simplified conditional estimating equations are compared with the conditional estimating equations with full and full covariance matrix .

Consider two-generation families composed of two parents and two children, . The proband is randomly selected from the four family members and is indexed by . A Weibull distribution is adopted for the onset time for all family members; , where is a binary variable with , , and we assume that , ; . Let , , and choose to give a median age of 45 years for disease onset for group with . The clinic entry time for the proband is normally distributed with mean 50 and variance 20, and families are recruited into the study only if their probands satisfy the selection condition . For non-proband j in the selected family i, let be the random age of contact, following for individuals in the first generation and for the individuals in the second generation, ; the age at contact for individuals in both generation are truncated at 90 years. We consider a Gaussian copula to induce an exchangeable within-family association for simplicity here, and let Kendall’s vary from 0 to 0.5 to reflect independence to strong within-family association. The second-order model with a Fisher transformation link function is simply , . The asymptotic variances of estimators based on conditional estimating equations in (2) are approximated by Monte Carlo simulation based on 20,000 samples.

Figure 1 shows the trends of asymptotic relative efficiencies of estimators under different conditional estimating equations as a function of the within-family association. It is apparent that the conditional estimating equations with full and full (G–W) lead to the most efficient estimators, and the efficiency gain is most appreciable for the association parameter. With the WPI matrix , adopting yields more efficient estimators than that using , especially when the within-family association is strong. This makes sense as the former utilizes additional information about from the conditional mean . The conditional estimating equations with the full and WPI matrix (G-WPI) perform worse than other approaches when the association is less than 0.45; see Fig. 1. This indicates that with a working covariance matrix, using the full derivative matrix increases the complexity, but does not improve efficiency; on the contrary, it leads to less efficient estimators. This is similar to the findings reported by Balemi and Lee [1] where they compare the performance of GEE1 and GEE2 estimators for clustered binary data.

Fig. 1

Asymptotic relative efficiencies of estimators under conditional estimating equations with full and WPI (), simplified G and WPI (), simplified G and WPI () compared with that under conditional estimating equation with full and full covariance matrix ; within-family dependence of disease onset times is induced by a Gaussian copula with exchangeable structure with Kendall’s varying from 0 to 0.5; , ,

Finite Sample Study of the Conditional Estimating Equations

Here we conduct a simulation study to assess the validity and finite sample performance of these four conditional estimating equations for family data from response-dependent sampling. The parameter settings are the same as in Sect. 3.1 and we let Kendall’s , and 0.4 for an exchangeable Gaussian copula. We also consider a more general case where the within-family association is induced by a Gaussian copula with structured correlation matrix. For the two-generation families composed of two parents and two children, let Kendall’s for parents, Kendall’s for parent–child, and Kendall’s for siblings to reflect the existence of both environmental and genetic effects on the age of onset. One thousand datasets of and 1000 ascertained families are generated, the four proposed conditional estimating equations are used for analysis, and the empirical properties of estimates of and are summarized in Table 1 for the exchangeable Gaussian copula and and for the structured Gaussian copula in Table 2. The performance of the estimators for the parameters of the baseline hazard was excellent under the (correct) Weibull specification in all settings and so we do not tabulate these results.

Table 1

Empirical properties of estimates under conditional estimating equations for family data from response-dependent sampling, where within-family association is induced by Gaussian copula with exchangeable structure; parametric margins and piecewise constant baseline hazards (3 pieces) are considered; ,

		Weibull margin								Piecewise constant margin^a
		β				γ_0				β				γ_0
τ	EE^b	EBIAS	ESE	ASE	ECP%	EBIAS	ESE	ASE	ECP%	EBIAS	ESE	ASE	ECP%	EBIAS	ESE	ASE	ECP%	NC %
m=200
0.0	G–W	0.002	0.113	0.115	95.5	-0.002	0.067	0.065	94.1	0.002	0.113	0.115	95.7	-0.000	0.070	0.067	93.4	0.2
	G-WPI	0.001	0.124	0.125	95.7	0.000	0.100	0.096	93.9	0.002	0.125	0.125	95.7	0.004	0.111	0.102	94.0	1.3
	G^I-WPI	0.001	0.113	0.115	95.7	-0.001	0.080	0.077	93.2	0.002	0.114	0.115	95.6	0.000	0.084	0.080	93.7	0.0
	G^II-WPI	0.001	0.113	0.115	95.7	-0.001	0.078	0.074	93.2	0.001	0.113	0.115	95.7	0.001	0.081	0.076	93.2	0.2
0.2	G–W	-0.006	0.109	0.109	94.5	0.002	0.106	0.105	94.6	-0.006	0.109	0.108	94.3	0.021	0.115	0.112	95.2	0.1
	G-WPI	-0.010	0.120	0.119	95.1	0.009	0.136	0.133	95.0	-0.011	0.120	0.118	95.1	0.050	0.163	0.151	93.8	2.2
	G^I-WPI	-0.006	0.109	0.109	94.2	0.005	0.118	0.117	94.8	-0.007	0.109	0.108	94.2	0.032	0.132	0.126	95.0	0.0
	G^II-WPI	-0.007	0.109	0.109	94.3	0.006	0.119	0.119	95.4	-0.007	0.108	0.108	94.1	0.031	0.133	0.128	95.8	0.1
0.4	G-W	0.001	0.092	0.094	96.0	0.008	0.152	0.154	95.6	-0.003	0.090	0.092	96.3	0.065	0.174	0.164	94.3	0.1
	G-WPI	0.001	0.103	0.104	96.4	0.015	0.173	0.175	95.7	-0.005	0.101	0.101	95.8	0.100	0.208	0.197	93.3	1.6
	G^I-WPI	0.001	0.093	0.094	96.0	0.013	0.170	0.173	96.6	-0.004	0.091	0.092	95.6	0.083	0.195	0.189	96.6	0.0
	G^II-WPI	0.000	0.093	0.094	96.1	0.016	0.183	0.186	96.0	-0.004	0.091	0.092	95.7	0.091	0.217	0.211	98.4	2.0
m=1000
0.0	G–W	-0.001	0.051	0.052	95.6	0.001	0.030	0.029	94.7	-0.002	0.049	0.052	95.9	0.001	0.030	0.030	95.3	0.0
	G-WPI	-0.001	0.055	0.056	95.4	-0.000	0.042	0.043	95.1	-0.002	0.054	0.056	95.6	0.001	0.046	0.046	94.9	0.0
	G^I-WPI	-0.001	0.051	0.052	95.8	-0.000	0.035	0.035	94.0	-0.002	0.049	0.052	96.0	0.001	0.036	0.036	94.9	0.0
	G^II-WPI	-0.001	0.051	0.052	95.8	0.000	0.034	0.034	93.9	-0.002	0.049	0.052	96.1	0.002	0.035	0.034	95.1	0.0
0.2	G–W	-0.004	0.050	0.049	93.2	0.002	0.048	0.047	95.2	-0.005	0.050	0.048	93.4	0.021	0.051	0.050	93.5	0.0
	G-WPI	-0.004	0.054	0.054	93.6	0.004	0.059	0.059	95.6	-0.005	0.054	0.053	93.5	0.039	0.068	0.067	92.7	0.0
	G^I-WPI	-0.004	0.050	0.049	93.4	0.003	0.053	0.052	95.0	-0.005	0.050	0.048	93.4	0.031	0.057	0.056	93.2	0.0
	G^II-WPI	-0.004	0.050	0.049	93.4	0.004	0.054	0.053	94.7	-0.005	0.050	0.048	93.4	0.030	0.058	0.057	92.8	0.0
0.4	G–W	0.000	0.042	0.042	95.0	0.003	0.068	0.068	95.0	-0.004	0.041	0.041	94.7	0.058	0.076	0.073	88.9	0.0
	G-WPI	-0.000	0.045	0.047	95.5	0.006	0.078	0.077	95.1	-0.007	0.044	0.045	95.3	0.088	0.092	0.088	83.7	0.0
	G^I-WPI	0.000	0.042	0.042	95.0	0.005	0.075	0.075	95.9	-0.005	0.041	0.041	94.7	0.076	0.085	0.083	86.7	0.0
	G^II-WPI	-0.000	0.042	0.042	94.9	0.005	0.080	0.080	96.0	-0.005	0.041	0.041	94.8	0.079	0.094	0.090	88.6	0.0

 Empirical properties are summarized based on replicates leading to convergence of the piecewise constant model. The percentages of replicates failing to converge (out of ) are indicated in the last column (all replicates converged for the approach based on the correct Weibull marginal model)

 G corresponds to with ; G corresponds to with also; WPI represents the working partial independence assumption with ,

Table 2

Empirical properties of estimates under conditional estimating equations for family data from response-dependent sampling, where within-family association is induced by Gaussian copula with structure, and ; parametric margins and piecewise constant baseline hazards (3 pieces) are considered; ,

		β				γ_0				γ_1				γ_2
Margins	EE^b	BIAS	ESE	ASE	ECP%	BIAS	ESE	ASE	ECP%	BIAS	ESE	ASE	ECP%	BIAS	ESE	ASE	ECP%	NC %
m=200
Weibull	G-W	0.001	0.103	0.104	96.3	-0.002	0.172	0.169	94.4	0.015	0.194	0.196	95.2	0.004	0.143	0.138	94.7	0.0
	G-WPI	0.001	0.115	0.114	95.3	0.006	0.239	0.233	94.2	0.015	0.267	0.262	94.3	0.001	0.196	0.191	94.2	0.0
	G^I-WPI	0.000	0.103	0.104	96.3	-0.001	0.216	0.210	93.3	0.018	0.262	0.256	94.8	0.004	0.193	0.187	94.6	0.0
PWC-3^a	G-W	-0.000	0.102	0.103	96.2	0.011	0.179	0.175	94.4	0.031	0.207	0.208	95.0	0.011	0.147	0.142	94.3	0.2
	G-WPI	-0.001	0.114	0.112	95.1	0.032	0.255	0.246	93.7	0.058	0.317	0.303	94.0	0.016	0.211	0.205	94.3	3.9
	G^I-WPI	-0.001	0.102	0.103	96.3	0.013	0.227	0.219	93.3	0.051	0.298	0.286	94.8	0.016	0.205	0.197	94.7	0.4
m=1000
Weibull	G-W	-0.002	0.045	0.047	95.9	-0.000	0.073	0.075	97.2	0.006	0.084	0.088	96.2	0.001	0.060	0.062	95.2	0.0
	G-WPI	-0.001	0.050	0.051	95.6	0.005	0.102	0.104	95.0	0.005	0.114	0.117	95.9	-0.001	0.084	0.085	95.8	0.0
	G^I-WPI	-0.002	0.045	0.047	96.0	0.002	0.091	0.094	95.6	0.007	0.111	0.115	96.0	-0.000	0.082	0.084	96.4	0.0
PWC-3^a	G-W	-0.003	0.045	0.046	95.6	0.014	0.075	0.078	96.5	0.021	0.092	0.093	95.6	0.007	0.063	0.063	95.0	0.0
	G-WPI	-0.003	0.050	0.051	95.3	0.030	0.108	0.110	94.4	0.040	0.135	0.135	94.2	0.011	0.090	0.092	95.9	0.0
	G^I-WPI	-0.004	0.045	0.046	95.7	0.019	0.095	0.098	95.0	0.036	0.126	0.128	94.9	0.010	0.087	0.088	95.5	0.0

Empirical properties are summarised based on replicates leading to convergence for the piecewise constant model. The percentages of replicates failing to converge (out of ) are indicated in the last column (all replicates converged for the approach based on the correct Weibull marginal model)

G corresponds to with ; G corresponds to with also; WPI represents the working partial independence assumption with ,

The results under the exchangeable Gaussian copula in Table 1 show that when the Weibull model is specified for the onset time distribution, the empirical biases are negligible for all conditional estimating equations; there is very slight finite sample bias for the association parameter when and the within-family association is strong (Kendall’s ). The empirical standard errors (ESE) agree with the average standard errors (ASE) based on the robust variance form, and the empirical coverage probabilities (ECP) of nominal 95% confidence intervals are in general within the acceptable range. Consistent with the theoretical results of Sect. 3.1, the greatest efficiency came from the conditional estimating equations with the full derivative matrix and full covariance matrix (G–W), followed by those with G and WPI matrix (G-WPI). The empirical performance of the conditional estimating equations with the full and WPI matrix is worse than the others, again in alignment with the conclusion based on Fig. 1.

When considering a more flexible marginal model with a piecewise constant (3 pieces) baseline hazard function, we set the cut-points at and 40. For the large sample size , performance was excellent for inference regarding and very good for the dependence parameter when the association was small; properties of the estimator of became worse with stronger within-family dependence, possibly as a result of the crude approximation of the piecewise constant hazard. While one might expect superior performance if more pieces were accommodated, convergence problems arose even with just three pieces under the smaller sample sizes for some replicates (typically less than 2.5%); the percentages of replicates failing to converge are reported in the last column and where necessary the properties of estimators from converged replicates are given. The G-WPI estimating equation always resulted in convergence. The convergence issues likely arose due to the right-truncated nature of the proband onset time and the severe censoring from a current status observation scheme of non-probands; these combine to yield little information to estimate the hazard function in small samples.

Under more general association structure, results under G-WPI estimating equation are not summarized because of high non-convergence percentage for such more general association structure. For the other three conditional estimating equations, their performance was again excellent under the correct Weibull model and again 100% of the replicates lead to convergence for and ; see Table 2. Empirical biases were generally small, there was a good agreement between the empirical and average robust standard errors, and the empirical coverage probability was generally within the acceptable range. Under the piecewise constant model, convergence rate was 100% when and the empirical properties of the estimators for and were good in such settings. When , performance remained good but with small finite sample bias and good empirical coverage probability.

Impact of Misspecifying the Dependence Structure

Limiting Bias Under Misspecified Conditional Estimating Equations

While standard GEE1 only requires correct specification of the marginal mean for consistent estimation of the marginal parameters, the conditional estimating equations require correct specification of the marginal distribution and the dependence structure for consistent estimation, even for the simplified conditional estimating equations G-WPI and G-WPI. As is often the case, the efficiency gains coming from the use of higher-order moments in the conditional estimating equations such as G–W come at the cost of poorer robustness. We explore the limiting behavior of estimators from misspecified models here based on large sample theory [31]. Specifically, if is an estimating function for based on a misspecified model, then the solution for asymptotically follows:as , where , , and is the solution to , where denotes an expectation taken with respect to the true distribution indexed by . Note that can be written aswhere and are the conditional expectations of and given under the true model. The expectation on the right-hand side of (6) is taken with respect to the remaining random variables . Of course, when the model is correctly specified, then but this is not the case more generally; we investigate the limiting bias of estimators under the misspecified model by examining .

Here we consider two-generation families composed of two parents and two children, and the proband is randomly selected from the four family members. The probands are recruited into the registry only if . We adopt the same parameter settings as in Sect. 3.1, but assume here that the true within-family association structure is induced by the Clayton copulawhere Kendall’s . The adopted estimating functions are misspecified in that the dependence structure is modeled based on a Gaussian copula with an exchangeable association structure. We consider the values of Kendall’s ranging from 0 to 0.5 to reflect independence to strong within-family dependence. We evaluate the limiting relative biases of estimators using Monte Carlo methods to take the expectation in (6) and solving the resulting equation.

From Fig. 2, we see that the conditional estimating equations with the full and WPI matrix is the most sensitive to misspecification. Although one might anticipate that the full and full covariance matrix (G–W) would be less robust than G-WPI or G-WPI, the asymptotic relative biases of estimators defined through G–W are in general no larger than those under G-WPI and G-WPI when Kendall’s is less than 0.3; the sensitivity of estimators from G–W to misspecification becomes more apparent, compared to those based on G-WPI and G-WPI, when Kendall’s is larger (i.e., ); G-WPI is slightly more sensitive to this form of misspecification than G-WPI. Furthermore, the asymptotic relative biases for under the conditional estimating equations are all relatively modest when Kendall’s is small to modest. If one is primarily interested in the estimation of , then the proposed conditional estimating equations are reasonably robust to misspecification of the copula function for modest Kendall’s , but the asymptotic biases of the dependence parameters are appreciable under misspecification of the dependence structure. This conclusion is analogous to those made regarding misspecification of the random effect distribution with response-dependent sampling [11, 13, 23]. We also conducted supplementary simulation studies demonstrating a good agreement between the finite sample and asymptotic biases in studies with 200 and 1000 families (see Supplementary Material).

Fig. 2

Asymptotic relative biases of estimators under conditional estimating equations when a Gaussian copula with an exchangeable structure is adopted for within-family dependence modeling; the true within-family dependence structure is induced by a Clayton copula;

Power Implications of Dependence Structure Misspecification

We next investigate the effect of dependence structure misspecification on the power of tests regarding covariate effects. Based on our previous findings regarding asymptotic relative efficiency and robustness, here we focus our attention on the preferred estimating functions G–W and G-WPI. We consider a test of versus , and let be the clinically important effect. When both the marginal and association models are correctly specified from (2), we haveas , where is the diagonal element in the robust covariance matrix corresponding to . Under a two-sided Wald test with significance level 100%, the required number of families to ensure 100()% power to detect is the smallest m satisfyingwhere and are the square roots of asymptotic variances of under the null and alternative hypotheses; and . is the percentile of standard normal distribution.

When the dependence structure is misspecified, the limiting value of estimators under the conditional estimating equations is (Sect. 4.1). Then based on (6), we can calculate the limiting values of under the null and alternative hypotheses when the dependence structure is misspecified, and denote them as and , respectively. Furthermore, we can show that under the null hypothesis the estimator based on the misspecified conditional estimating equations satisfiesas , and under the alternative hypothesiswhereHence, the asymptotic properties of can be determined by considering the corresponding component of . When the copula model is misspecified, the actual power of such two-sided Wald test of versus , at the clinically important effect given sample size m and significance level , iswhere and are the square roots of the diagonal elements of and , respectively, corresponding to .

Here we report on an asymptotic study to examine the effect of copula misspecification on the power. Assume that each family consists of two parents and two children, and the proband is randomly selected from the four family members. As before we presume that the families are recruited to the study only if . The parameter settings are the same as in Sect. 3.1 but we consider two specific scenarios: (i) the true within-family association structure is based on a Gaussian copula with an exchangeable association structure and (ii) the true within-family association structure is based on a Clayton copula (7); in both cases, we set Kendall’s . At the design stage, we adopt a Gaussian copula with an exchangeable association structure for the within-family dependence, and let Kendall’s . We therefore only consider the case in which the form of the dependence structure is misspecified. In this setting, we calculate the required sample size to achieve 80% power to reject at by (9), where is obtained from the Gaussian copula. The minimum numbers of families are 420 and 422 based on estimating equation G–W and G-WPI, respectively. Under these sample sizes, the actual power of such a design can be computed by (12) for the values of ranging from 0 to . The power curves are plotted in Fig. 3 from which we infer that when the association model is correctly specified, tests based on the conditional estimating equations G–W and G-WPI have the desired power at the clinically important effect; as expected, the power decreases when the true value of approaches 0. When the copula is misspecified (i.e., the true dependence structure is set by a Clayton copula but a Gaussian copula is used for sample size calculation), tests based on both conditional estimating equations lead to a loss in power, with a greater loss in power under G–W compared to G-WPI. This is reasonable since the G–W estimating equations exploit information from higher-order dependencies more than G-WPI, which is less robust than the latter. In summary, based on the comprehensive investigation of these conditional estimating equations in terms of efficiency and robustness, estimating equation G-WPI is suitable in the absence of information about the association structure, but if information is available about the structure, estimating equation G–W could be adopted to achieve higher efficiency.

Fig. 3

Power curves of a two-sided Wald test for under conditional estimating equations G–W and G-WPI when the within-family dependence structure is correctly specified or misspecified; true within-family dependence is induced by Gaussian copula with exchangeable structure or Clayton copula, and adopted family dependence structure in the design stage is Gaussian copula with exchangeable association; Kendall’s ,

Application to The Psoriatic Arthritis Family Study

Hereditary factors are thought to be important in psoriatic arthritis, as some studies have suggested that close blood relatives of affected individuals are at higher risk of developing the disease compared to the general population. Interest therefore lies in characterizing the effect of genetic markers on the risk of disease; we consider four human leukocyte antigen (HLA) markers reported in the literature as being associated with psoriasis or psoriatic arthritis including HLA-B8, HLA-B27, HLA-C6, and HLA-C12. Characterizing the nature of the within-family association structure can also provide useful insight into the genetic basis of the disease. Particular interest lies in assessing the “parent-of-origin” effect; preliminary evidence suggests that there may be a stronger risk of paternal transmission, over maternal transmission, of risk of disease; we refer readers to Pollock et al. [24] for associated results based on binary analyses.

Here we consider an application to the motivating family study on the genetic basis of psoriatic arthritis conducted in the Centre for Prognosis Studies in the Rheumatic Diseases at the University of Toronto. One hundred and sixty-nine families composed of 2 to 7 members, including the proband, were recruited for study. The date of disease onset is available for probands from the clinic registry but only the disease status of other individuals is available when they are examined, yielding current status data. A Weibull model is adopted for the marginal distribution of the PsA onset time with survivor function where , , and . A flexible model for the within-family dependence is formulated based on a Gaussian copula with different pairwise dependencies between parents (), between siblings (), between a father and his child (), and between a mother and her child (). This can be formulated in terms of a second-order regression model given bywhere , , and . The hypotheses and are the basis of a test regarding the parent-of-origin question. There are only 8 pairs of parents, which leads to insufficient data to estimate the intercept in (13) and so we constrain that parameter to be zero with the implicit assumption that there are no environmental determinants of PsA.

Table 3 summarizes the results with the top half obtained from the full derivative and covariance matrices (G–W) and the bottom half reporting the results from G-WPI. A model with no HLA covariate is given in the first column followed by four univariate models, with the last column containing results from a multivariate model including all four markers. The estimates for the association parameters are given in terms of and the three Kendall’s parameters.

Table 3

Estimates of analyses of HLA markers and time to disease onset based on conditional estimating equations using response-biased psoriatic arthritis family data; associated standard errors are in the parentheses

	Univariate models					Multivariate model
	No HLA	HLA-B8	HLA-B27	HLA-C6	HLA-C12
Full G and Full W (G-W)
logλ	-5.461 (0.531)	-5.353 (0.455)	-5.573 (0.609)	-5.650 (0.590)	-5.342 (0.464)	-5.544 (0.543)
logκ	1.347 (0.086)	1.347 (0.086)	1.348 (0.086)	1.349 (0.086)	1.346 (0.087)	1.349 (0.086)
β_B8	-	-0.483 (0.709)	-	-	-	-0.195 (0.751)
β_B27	-	-	1.161 (0.755)	-	-	1.308 (0.771)
β_C6	-	-	-	0.397 (0.593)	-	0.458 (0.625)
β_C12	-	-	-	-	0.788 (0.769)	0.982 (0.830)
γ_1	0.700 (0.242)	0.660 (0.224)	0.741 (0.261)	0.777 (0.242)	0.620 (0.234)	0.695 (0.244)
γ_2	0.458 (0.264)	0.414 (0.242)	0.470 (0.281)	0.536 (0.259)	0.378 (0.264)	0.423 (0.269)
γ_3	0.261 (0.291)	0.220 (0.268)	0.300 (0.310)	0.343 (0.285)	0.173 (0.289)	0.252 (0.287)
τ_ss	0.337 (0.107)	0.318 (0.101)	0.354 (0.114)	0.370 (0.104)	0.301 (0.106)	0.334 (0.108)
τ_fc	0.225 (0.125)	0.204 (0.116)	0.231 (0.133)	0.262 (0.121)	0.187 (0.127)	0.209 (0.129)
τ_mc	0.130 (0.143)	0.110 (0.132)	0.149 (0.151)	0.170 (0.138)	0.086 (0.143)	0.125 (0.141)
G^I and WPI W (G^I -WPI)
logλ	-5.037 (0.326)	-5.007 (0.313)	-5.105 (0.358)	-5.034 (0.325)	-5.052 (0.317)	-5.105 (0.331)
logκ	1.337 (0.087)	1.338 (0.087)	1.339 (0.087)	1.337 (0.087)	1.338 (0.088)	1.340 (0.087)
β_B8	-	-0.490 (0.644)	-	-	-	-0.325 (0.669)
β_B27	-	-	0.849 (0.587)	-	-	0.969 (0.580)
β_C6	-	-	-	-0.037 (0.492)	-	0.060 (0.537)
β_C12	-	-	-	-	0.767 (0.644)	0.911 (0.667)
γ_1	0.448 (0.232)	0.447 (0.234)	0.481 (0.237)	0.447 (0.238)	0.437 (0.227)	0.465 (0.238)
γ_2	0.208 (0.237)	0.199 (0.234)	0.208 (0.238)	0.206 (0.241)	0.195 (0.241)	0.180 (0.244)
γ_3	-0.036 (0.249)	-0.039 (0.246)	-0.005 (0.250)	-0.037 (0.256)	-0.059 (0.251)	-0.036 (0.254)
τ_ss	0.220 (0.110)	0.220 (0.111)	0.236 (0.112)	0.220 (0.113)	0.215 (0.108)	0.228 (0.113)
τ_fc	0.104 (0.117)	0.099 (0.116)	0.104 (0.118)	0.103 (0.119)	0.097 (0.119)	0.090 (0.121)
τ_mc	-0.018 (0.124)	-0.019 (0.123)	-0.002 (0.125)	-0.019 (0.128)	-0.030 (0.125)	-0.018 (0.127)

Based on the model with no HLA covariates, we find (95% CI: 0.113, 0.528; p value = 0.002), indicating highly significant association between siblings in the disease onset time. The father–child association is lower at (95% CI: −0.030, 0.452) and not quite statistically significant (p value = 0.072). For the mother–child association, we find (95% CI: −0.153, 0.393) which is weaker still and insignificant (p value = 0.364). A test of the parent-of-origin hypothesis based on yields a Wald statistic of 1.435 (p value = 0.151). As this is not statistically significant at the 5% significance level, there is insufficient evidence to claim a statistically significant “parent- of-origin” effect. The results are broadly comparable for the HLA regression analyses based on the other conditional estimating equation (G-WPI). For the association parameters, the estimates are somewhat lower with (95% CI: −0.003, 0.423; p value = 0.046), (95% CI: −0.128, 0.324; p value = 0.378), and (95% CI: −0.256, 0.222; p value = 0.886). The Wald statistic of 1.682 (p value = 0.092) does not suggest a “parent- of-origin” effect.

The large sample theory we develop can be used to plan a future family study and it is possible to calculate how many families would be required to ensure adequate power to test the parent-of-origin hypothesis in a future study. In a new study, we may consider recruitment of families of members of the registry and presume that the distribution of family members, ages at assessment, and other factors are similar in the new study. We use the sample size formula similar to (9) but for and determine that 627 families would be required to ensure 80% power to detect a significant difference between the father–child and mother–child associations using estimating function G–W when the true effects correspond to those seen in the first column of Table 3. The current study therefore appears to be grossly under-powered to formally test the parent-of-origin hypothesis.

None of the HLA markers were shown to have a significant effect on the time to the onset of PsA. Based on the G–W estimating equations, there is a trend toward a reduction in risk with HLA-B8 and a trend toward an increased risk with the presence of each of the other HLA markers.

Discussion

Estimating functions have been developed to model the nature and extent of within-family dependence in disease onset times from family studies under response-dependent sampling. A novel aspect of this work is the formulation of the dependence measures on the basis of the disease onset time and the recognition that the available data on family members are handled more naturally as current status data rather than binary data. This approach utilizes all available data from probands and their relatives in assessing the association between age of onset and covariates and in evaluating the association structure of age of onset among family members. The biased sampling scheme typically employed in family studies is addressed by the use of conditional estimating equations where the conditioning event reflects the selection criteria. Several specific estimating functions within the class proposed are assessed in terms of efficiency and robustness; these results complement the standard results of second-order estimating functions since all moments in the proposed equations are conditional. We also outline how sample size requirements for family studies can be assessed based on this framework to ensure that power objectives are met. Code for solving the conditional second-order estimating equations (1) and for obtaining the variance estimates of Sect. 2.2 are available at Github https://github.com/Yujie-Zhong/CGEE2.

We have focused on the use of estimating functions for the analysis of family data in part because the likelihood can be challenging to compute when the size of the family is large. Nevertheless, some assessment of the loss of efficiency in comparison to this optimal approach would be worthwhile. The validity of the proposed conditional second-order estimating equations hinges on correct specification of the dependence structure, a requirement that is analogous to the need for correct specification of the mixing distribution in random effects models for data obtained based on a response-dependent sampling scheme [23]. Assessing model adequacy is best done by testing for the need for model expansion; this could be carried out by testing the need for more cut-points in the baseline hazard function to accommodate a more flexible hazard function, or the need to test for a more general dependence structure. In the present setting, the dependence structure is most easily formulated by selecting a working copula model for the joint distribution of the onset times in the population. If this dependence structure is misspecified, inconsistent estimates are obtained, and we examine their consequences in Sect. 4 to make recommendations on the use of a particular derivative and working covariance matrix. The properties of estimators under model misspecification can be explored using large sample theory [31], but these will be influenced by response-dependent sampling schemes and so a more general study of the effect of misspecification in this framework represents an important area for further research.

We have restricted our attention to parametric models for the onset time distribution. Natural extensions would be to introduce non-parametric or semi-parametric methods for estimating the marginal distributions. In the latter case, one can look at multiplicative Cox models, accelerated failure time models, and Aalen’s additive model, among many other methods. Joint estimation based on the most general conditional estimating equation can be challenging in this setting, but two-stage estimation procedures may be feasible; this is an area of current research. The preliminary work based on the piecewise constant baseline hazard model, however, suggests that studies may need to recruit a lot of families if the incidence rate is low to estimate the marginal onset time distribution. If the disease onset times are available for all or even some of the non-probands found to have the disease, these data could help in estimation; the estimating equations we present can be modified in this case to incorporate such data. Auxiliary samples can also be useful to enhance inferences.

While there is an increasing amount of attention given to the use of disease onset time as a basis for modeling within-family dependence, there remain challenging issues that warrant further attention. The primary challenge is in quantifying dependence in the presence of the competing risk of death [12, 26]. The classical illness–death process is a natural framework for modeling the occurrence of disease in individuals who are at risk, and generalization of this set-up to model within-family dependence is an area warranting attention. This issue is not unique to analyses based on disease onset times; when current status data are treated as binary data, the requirement that individuals are alive at the time of contact is ignored.

Below is the link to the electronic supplementary material.

Supplementary material 1 (pdf 53 KB)