Optimal Design Strategies for Sibling Studies with Binary Exposures

1 Introduction

In recent years, epidemiological studies involving siblings have become increasingly common [1–7], and this design has now also been recognized in popular textbooks [8]. Sibling designs can provide better control over confounding by unmeasured family-level risk factors shared by the siblings (e.g. genetic, environmental, or socioeconomic) compared to other sampling strategies [9]. In addition, they can enhance study efficiency by reducing extraneous variability. Sibling studies are sometimes limited to two individuals per family, as in the case of twin studies [10], but more commonly comprise sibships of varying sizes, including singletons [11].

To estimate exposure effects with correlated data, two popular approaches are generalized linear models (or marginal models) using generalized estimating equations (GEE) to handle within-family correlations and random effects models (or mixed effects models). Marginal models have a population-average interpretation, whereas mixed effects models capture variation across families and allow for a family-specific interpretation. There has been extensive discussion in the literature of the relative merits of marginal models compared with mixed effects models [12, 13]. A recent study by Hubbard et al. [14] advises that mixed effects models should be approached with caution as they rely on unverifiable distributional assumptions. Specifically for sibling studies, a number of approaches have been employed to date including between–within (BW) models with fixed or random intercepts [15], conditional logistic regression models [16] for binary outcomes, and marginal models [9, 17]. Several review articles [9, 15, 18] provided detailed comparisons of the various approaches. In the context of fixed effect BW models and marginal models, Sjölander et al. [18] and Frisell [9] pointed out that the estimation and interpretation of the model parameters depends on the shared and non-shared confounders of sibling sets. Specifically, within-sibling estimates from BW models have larger bias in the presence of non-shared confounders relative to those from marginal models [9]. Furthermore, within-sibling estimates will have larger bias if the within-family correlation of exposure is higher than that of confounders.

Still missing from the literature, however, is a comprehensive treatment of optimal design strategies for sibling studies and easy access to the calculation of study power in specific settings. Ideally, study power should be readily available for any combination of exposed and unexposed individuals in sibships of varying sizes, and take within-family correlations into account, for either continuous or discrete study outcomes. To address this need, we here present design strategies to optimize study power based on marginal models using GEE to handle the within-family correlations. GEE provides an explicit way of handling within-family correlation, which for the purpose of power calculation is especially convenient. As far as we know, explicit power calculation methods for mixed effects models in the case of binary outcomes have only been developed under fixed alternatives [19], or rely on simulation [20] and the use of unverifiable distributional assumptions [14]. In the GEE setting we are able to explicitly calculate the asymptotic power (based on results in Li and McKeague [21]) and avoid the use of simulation methods. Throughout this paper, we use E and C to denote exposed and control (where “control” refers to unexposed) status, respectively. So an EC sibling pair (or an EC family) denotes two siblings with opposite exposure status, and an EE (or CC) sibling pair denotes two exposed (or unexposed) siblings.

The paper is organized as follows. Preliminary materials setting up the statistical model appear in Section 2. The proposed design strategies are developed in Sections 3 and 4. Section 5 contains discussion.

2 Statistical model

In this section, we formulate the statistical model that is used throughout the paper. Let yij be the outcome, xij the exposure of interest for the jth sibling in the ith family, and μij=E(yij|xij) be the conditional mean of yij given xij. A marginal model [9, 12] to estimate the effect of exposure can be written as

where g(.) is logit link for binary outcome and identity link for continuous outcomes, parameter λ is the intercept and parameter β measures the association between exposure and outcome. The parameters associated with this model are estimated by the GEE approach [21, 22] which requires specifying a working correlation matrix. We assume that siblings have positive correlations in terms of the outcome of interest.

We are interested in testing the hypothesis that a specific exposure is associated with a specific outcome of interest, i.e. H0:β=0 vs H1:β≠0. Our results are based on the use of a quasi-score test (rather than a Wald test) [21]. Quasi-score tests under marginal models have been studied extensively [23, 24] and provide a sound alternative to the score test when a score function is not available (the likelihood function is often intractable for correlated data). Unlike a score statistic, a quasi-score test is constructed from an estimating equation and an associated sandwich-type variance estimator, which results in a test that is asymptotically equivalent to a quasi-likelihood ratio test. The Wald test is a popular test, but Hauck and Donner [25] reported that the Wald statistic in logistic regression decreases under the alternative with increasing effect size because of increasing variance of the estimator of the regression parameter. This “aberrant” behavior of the Wald statistic raises difficulties for the study of optimal design strategies and is beyond the scope of the present paper.

Let m be the number of families in the study, α the type I error rate, and x1,1−α2 the 100(1−α)th percentile of the chi-square distribution with 1 degree of freedom. According to Li and McKeague [21], the test statistic asymptotically follows a chi-square distribution with 1 degree of freedom and the power of detecting a specific β=β1 is the probability P[χ12(νm)≥χ1,1−α2], where χ12(νm) is a non-central chi-squared random variable with 1 degree of freedom and the non-centrality parameter νm is given in eq. (2) in the Appendix. It is straightforward to see that the statistical power increases with the non-centrality parameter. With this formula, we can calculate study power for family studies with varying sibship sizes and allocation schemes. As study power is uniquely determined by the non-centrality parameter, a comparison of non-centrality parameters will provide the differences in statistical power for alternative design strategies. In the following sections, various common design issues that arise in sibling studies will be explored to provide recommendations for optimizing power.

3 Optimal design strategies for studies with a continuous outcome

We first evaluate optimal design strategies for studies with a continuous outcome for the following scenarios: (1) studies including families with one singleton or two siblings; (2) studies adding siblings to existing families; and (3) studies in which all families have the same number of siblings. The within-family correlation of the outcome is denoted by ρ throughout this paper. R code for calculating statistical power and minimum detectable effect size is provided in Sections S3, S5, S8, and S10 of the Supplemental Materials.

3.1 Studies including families with one singleton or two siblings

We start with a simple situation where up to two subjects are available from a single family and identify exposed individuals in each family by the letter “E” and unexposed (control) individuals by the letter “C”. With this notation, the study may include the family structures EC, EE, CC, E, and C. A first question then is whether to recruit one or two individuals per family. As mentioned earlier, the idea is simply to compare the non-centrality parameter νm for the two competing strategies. Let mec,mee,mcc,me, and mc denote the number of EC, EE, CC, E, and C families in a study, respectively. For continuous outcomes, the link function g(.) in model (1) is the identity link. That is μij=λ+βxij, and the conditional variance is given by var(yij|xij)=σ2 which does not depend on the mean μij. Notice that the conditional variance σ2 is assumed to be constant across the exposure levels. The non-centrality parameter in this case is given in eq. (3) in the Appendix. This allows for the comparison of various design strategies with different combinations of EE, CC, EC, E, and C families. These comparisons (in section “Non-centrality parameter for continuous outcomes involving EE, CC, EC, E and C families” in the Appendix) show that for a given number of subjects, the non-centrality parameter νm will be maximized when only EC families are included for the study. If this ideal design strategy is fully implemented, the non-centrality parameter reduces to

νm=β12σ2mec2(1−ρ),

where β1 is the alternative value and 1−ρ is the “design effect” in this case. As νm is an increasing function of ρ, study power will increase as ρ increases; largest study power is therefore obtained when outcomes are most highly correlated within sibships.

Sometimes in practice, not every exposed subject will have an unexposed sibling, in which case the above “ideal design strategy” cannot be implemented. However, from comparisons of the non-centrality parameter, it can be seen that recruiting two E (or C) families will be more efficient than recruiting a single EE (or CC) family. This also makes sense intuitively given that the positively correlated EE (or CC) siblings include redundant information whereas E (or C) singletons are independent.

From these observations, the following recruiting principles can be distilled:

1. Recruit EC sibling pairs whenever possible,

2. When (a) is not possible, recruit an equal number of individuals from E and C families.

If these recruiting principles are followed, a mixture of E, C, and EC families will be recruited for study (unless the ideal recruiting strategy is possible and the study can be limited to EC families). The non-centrality parameter νm is then given by

νm=β12σ2mec(mec+me+mc)1−ρ2+mcme2mec1+ρ+me+mc.

Notice that νm is still an increasing function of the correlation ρ, associated with the benefit gain from EC sibling pairs. Figure 1 illustrates the ideal strategies in relation to other strategies in terms of statistical power for a specific situation with an example of 80 subjects.

Figure 1

Power of sibling (solid) and singleton (dashed) designs; n = 80 individuals with 40 exposed, 40 unexposed, as the number of exposure–control sib pairs increases from 0 to 40

3.2 Adding siblings to available families

In typical epidemiologic studies, the number of exposed subjects is limited, but the number of unexposed controls is not. Let us therefore consider the situation in which some EC, E, and C families have been recruited for the study, and it is feasible to recruit additional control subjects. We could then recruit sibling controls either (a) from E families to form EC families, or (b) from EC families to form ECC families, or (c) to add C families.

Following the principles outlined above, attempts should first be made to recruit additional unexposed controls from existing E families to form EC sibling pairs. This will give a better within-sibship exposure contrast than forming ECC families or adding C families.

However, if additional controls are not available for E families, then should the extra controls be recruited from existing EC families to form ECC families or should C families be added? To compare the latter two scenarios we compare the corresponding values of νm. The expression of the non-centrality parameter is given in eq. (4) in the Appendix. The comparison shows that recruiting additional controls to form ECC families will generate greater statistical power than recruiting new C families, as illustrated in Figure 2, where extra controls are added to 40 EC pairs.

Figure 2

Additional power from adding extra unexposed individuals to existing EC families (solid) and as singletons (dashed) starting from 40 exposure-control sib pairs

4 Optimal design strategies for studies with a binary outcome

It suffices to consider binary outcomes (e.g., a disease outcome) with a prevalence of less than 0.5 in both the exposed and unexposed groups because if the disease prevalence is greater than 0.5, “non disease” status can be used as the outcome with a prevalence less than 0.5 for mathematical equivalence. R code for deriving optimal design strategies and calculating power and minimum detectable relative risk is provided in Sections S1, S2, S4, S6, S7, and S9 of the Supplemental Materials.

4.2 Adding siblings to available families

We next consider for binary outcomes the scenario with a fixed number of exposed study subjects available from EC, C, and E families and the possibility to recruit additional unexposed control subjects. Again, our options will be to recruit sibling controls either (a) from EC families to form ECC families, or (b) from E families to form EC families, or (c) to add C families. From the principles above it follows that one should first try to recruit additional unexposed controls from existing E families to form EC families. This will provide a better contrast of exposure status compared with adding C families or forming ECC families. The explicit formula of the non-centrality parameter is given in eq. (6) in the Appendix. With this formula, different recruitment strategies can be compared. As shown in Figure 3, where extra controls are added to 200 EC pairs, recruiting sibling controls from EC families will provide more statistical power compared to recruiting additional C families.

Figure 3

Additional power from adding extra unexposed individuals to existing EC families (solid) and as singletons (dashed) starting from 200 exposure-control sib pairs

5 Discussion

In this paper, we have investigated the design efficiency of sibling studies with binary exposures using a novel GEE-based [21] approach to calculate study power and estimate required sample sizes. Correlations of (continuous or binary) study outcomes between siblings are taken into account. Our results are obtained by maximizing explicit expressions for the non-centrality parameters of the (chi-squared) limiting distributions of quasi-score test statistics under local alternatives. The optimal design strategies for continuous and binary study outcomes are found to differ while sharing some common elements.

For studies with at most two siblings per family and a continuous study outcome, the optimal design strategy is to recruit EC sibling pairs. When the same number (more than two) of siblings is to be recruited from all families, the optimal strategy is to balance the exposure status within each family as evenly as possible. However, for studies with a binary outcome, the variation of the outcome is different across the two exposure groups due to its dependence on the mean (or prevalence). This needs to be taken into account and results in a more complex form for the optimal design strategies. Unknown values of p0 and p1 will make implementation of the proposed optimal strategies challenging. In practice, the underlying prevalence will almost always be unknown and will have to be estimated from the literature. In the absence of any empirical estimates, design strategies could be evaluated for a plausible range of values and then conservative estimates selected. Besides the aforementioned differences, studies with continuous and binary outcomes also share some common elements. Recruiting an EE or a CC sibling pair is less efficient in terms of study power compared to recruiting two E singletons or two C singletons. Furthermore, adding an additional control to E or EC families to form EC or ECC families is more efficient than adding a new C family with only a singleton.

Our results show that all study subjects in family studies contribute to study power, albeit in different degrees in different settings. The efficiency of a family study will therefore be compromised by the exclusion of any study subjects, even if these subjects are not easily compared to siblings who differ in exposure status. As pointed out by Frisell [9], even results from sibling designs have to be interpreted with caution as they can be biased by confounders not shared by the siblings. Given the availability of analytic techniques that account for correlated and non-paired data, approaches that do not include all available study subjects [5, 6, 27, 28] should be avoided in the interests of validity and study power.

We have limited our focus to scenarios involving families with no more than three siblings, or the same number of siblings in each family; this should be sufficient for evaluating the most commonly used sibling designs. R code is included. Optimal recruiting strategies for settings with more complex family structures are not presented, because the maximization of the non-centrality parameter becomes highly computer intensive, and it would be beyond the scope of the article to present the results in any detail.

Our approach was developed and illustrated for cross-sectional sibling studies without repeated measurements. In a longitudinal sibling study with repeated measures, the key difference is that the correlation structure will be more complex because of the added within-subject correlation. The theoretical basis for the approach developed here [21] does not assume a specific correlation structure, so the formula provided in section “General formula for the non-centrality parameter” in the Appendix could potentially be used to provide guidelines for longitudinal sibling studies as well.

We have focused on the quasi-score test, but in practice the Wald test is popular, despite its limitations for binary outcomes mentioned earlier. The optimal strategies proposed in this paper also apply to the Wald test in the case of continuous outcomes. An interesting topic for future research would be to compare the quasi-score, quasi-likelihood, and Wald tests in terms of optimal design strategies for binary outcomes.

Our analysis defines study efficiency in terms of statistical power, ignoring any considerations relating to the relative cost of recruiting either exposed or unexposed siblings in different family settings. It can already be seen, however, that in common situations where study outcomes are positively correlated between siblings, the recruiting of EE or CC families might be easier and more cost effective compared to the recruiting of EC or E or C families, although this will be less efficient in terms of statistical power alone for a given sample size. When required, the trade-off between study cost and statistical efficiency could further be explored using the formulas in this paper given the costs of recruiting siblings and singletons.