In the epidemiological research, it is important that the collected data are translated into interpretable results which can be easily communicated to clinicians. The need for “translatable” evidence from research studies is of prime importance in the evaluation of clinical interventions, because they hold the potential to immediately influence the course of patient treatment. When evaluating these studies, the examination of “Effect Size” or (ES) can be a useful measure of the comparative efficacy of the treatment under investigation. In randomized clinical trials, an effect size estimate quantifies the direction and magnitude of an effect of an intervention.
When exposure and disease risk are measured on a binary scale, several measures of effect size are in current use  . The odds ratio (OR), the relative risk (RR), and the population attributable risk (AR) are the most commonly used measures of effect size in clinical as well as analytic epidemiology.
The concept of AR was introduced in  and is a widely used measure of the amount of disease that can be attributed to a specific risk factor. The AR combines the relative risk (RR) and the prevalence of exposure to measure the public health burden of a risk factor by estimating the proportion of cases of a disease that would not have occurred if we remove the risk factor.
The concept of AR and its statistical characteristics have been reviewed in  and in several publications   . Statistical inferences on AR require the availability of data from subjects randomly assigned to intervention groups. However, when the sampling strategy involves aggregate or clusters of individuals, adjusting for the effect of intracluster correlation is essential in order to conduct valid statistical inferences  and the references therein. However, the statistical pro- perties of estimators of AR when clusters are sampled have not yet been fully explored. The fundamental objective of our work is to fill the gap of performing statistical inference on AR under the clustered binary data situation.
In this paper, we obtain the variance of the estimated AR under cluster sampling, focusing on cohort and cross-sectional designs. In Section 2, we construct an AR estimator, and in Section 3, we derive its large sample variance adjusted for the intracluster correlation (ICC). In Section 4, we consider the split cluster design, and describe situations where we compare two correlated AR parameters. In Section 5, we conduct a Monte-Carlo experiment to evaluate the empirical power of Wald’s test on the null hypothesis of equality of two correlated attri- butable risk parameters. At the end of each section, we provide an example.
2. AR from Cluster Sampling
We start with a parallel group design where k clusters are exposed to a specified risk factor, and l clusters are not exposed, as in the data layout given in Table 1.
In Table 1, we assume that clusters have been selected at random from a well-defined population of exposed individuals, where the cluster has units. All individuals in this sample are assumed to be exposed to the risk factor. Let with and 0, denoting positive and ne- gative responses corresponding to the presence of the exposure with . Similarly we assume that clusters have been selected from the population of unexposed individuals, where the cluster has units. All units in the clusters can serve as controls assuming the
Table 1. Typical data layout for clustered data in two groups: exposed and unexposed.
absence of exposure. In the unexposed clusters, let with and 0 denote positive and negative responses with . Furthermore, let and denote respectively the total number of events in the exposed and non-exposed groups; provided that the misclassification error is zero. Therefore, conditional on , has binomial distribution with parameters . Similarly, conditional on has binomial distribution with parameters . To introduce a within cluster correlation, we assume that follows a beta distribution with probability density function (pdf) given in (1).
and that follows a similar beta distribution where pdf is denoted . The effect of the intracluster correlation among the responses may be accounted for as follows.
Under the transformations, , , the mean and variance of are given respectively by and
Similarly, and . Therefore, we have , . Consequently, the unconditional distribution of is beta binomial with, , and . Similarly; , and
It should be noted that the beta distribution assumptions imposed on the model parameters is not necessary, and one may adopt a quasi-likelihood set-up, by specifying the first two moments for in  and  . This set-up, would lead to the same expressions for and . The reason for introducing the beta distribution here is that while it serves as mechanism to create within cluster correlation, it will form the basis for generating data from beta binomial distribution using Monte-Carlo simulations in Section 5.
The parameters and are respectively interpreted as the within cluster correlations among all pairs of scores in the group of exposed and unexposed. We may obtain consistent estimators of and on replacing the parameters, , , and with appropriate estimators from the data as will be shown in the next sections. We shall now construct unbiased point estimators for the parameters and .
From  and  , we have has , . Similarly, has ,
, where , ,
, and . Clearly and are unbiased point estimators for respectively.
The data, under the above set up can then be summarized in a 2 ´ 2 table as shown in Table 2.
Formally, the AR is defined in  as:
where is the percentage of disease in the population, and is the percentage of disease in the population in the absence of exposure to the risk factor. Levin  defines the Attributable Risk as “the amount of disease that can be attributed to a specific risk factor”.
Using Bayes theorem, and from [  ; page 73], Equation (2) may be written as:
Here; RR is the relative risk or the risk ratio, and is the risk of expo-
sure. The RR is defined by . In terms of population parameters,
the AR as defined in (3) is equivalent to:
Under the transformation, , we get . We shall use this
transformation to facilitate the derivation of the large sample variance of AR.
The sample estimator of AR, is obtained using the data in a 2 ´ 2 cross classification as given in Table 2.
Epidemiologists use this statistic quite frequently to assess the consequences of an association between a binary outcome of interest and exposure to a risk factor . The total number of observations in the non-exposed and the exposed groups are given respectively by and , assumed fixed.
For a cross sectional or cohort study designs the estimator is from 
Table 2. Disease-exposure cross classification in a 2 ´ 2 table.
Following  , we shall derive the asymptotic variance of .
We first write, , and, ,
where , , , and .
Using the delta method  we can show to the first order of approximation that:
A consistent estimator of may be obtained on replacing the para-
meters by their moment estimators. An confidence interval on AR is thus given as:
The moment estimators of the intraclass correlations are obtained separately from the groups of exposed and unexposed clusters. The moment estimator of is given by:
Similar expressions for the (MSW, MSB) are obtained for the clusters of unexposed. The quantities (MSW, MSB) are estimated from the one-way ANOVA model when the responses are measured on the binary scale. For details the readers are referred to  .
We now consider two examples, the first is from data arising from a cross sectional study and the second example is on data from a randomized prospective trial.
Example 1: Cross-Sectional Study: The effect of consanguinity on congenital heart defects (CHD).
The Saudi Arabian CHD registry  was established in 1998, and by 2003 the registry evolved into Multi-Institutional research collaboration. The prime aim of this institution is to develop a registry whereby data from major referral hospitals across the country can provide patient information.
The participating hospitals are from regions that cover the country making the registry a nationwide data repository for the Kingdom of Saudi Arabia [Congenital Heart Disease Registry 2013]. The present example uses data on a major congenital heart disease; Patent Ductus Arteriosus (PDA). The incidence of PDA has been reported to be approximately, 1 in 2000 births, which accounts for 5% to 10% of all congenital heart diseases with female to male ratio of almost 2:1  . The PDA was found to occur with increased frequency in several genetic syndromes, with precise mechanisms resulting in persistent PDA not yet clear   .
Arab countries are notorious for consanguineous marriages, with first cousin types being the most common. For example in Jordan the prevalence of consanguinity was reported in  as 51.3%, Yemen, 40% as reported  , and almost 57% in Saudi Arabia as reported   . More recently, a survey of Saudi families conducted in  , estimated the prevalence of consanguinity to be as high as 56%.
For illustrative purposes of the methodologies presented in this section, we sampled two children from the registry whose mother are non-diabetic, with maternal age less than 40 years. Each sampled child was classified according to the presence/absence of PDA, and the type of parental consanguinity (exposure variable). Therefore, for the exposed (children from consanguineous marriages restricted to first degree cousin) and non-exposed (children from non-consangui- neous marriages) the cluster size is . The data are presented in Table 3.
Direct applications using Equations (5), (8), (9), and (10) we get:
, , ,
The square root of Equation (6) gives , and the 95% confidence interval of AR is: −0.104 < AR < 0.210.
The AR estimate is interpreted as follow: if among infants born with CHD, gi- ven that PDA among infants with CHD is a preventable event, then prohibiting first degree relatives’ marriages will reduce the chance of having PDA by 6%.
Example 2: Prospective Cohort study (Weil’s data)
The data in this example was given first in  taken from  , and gives the
Table 3. Consanguinity and PDA.
results from an experiment comparing two treatments. One group of 16 pregnant female rats was fed a control diet during pregnancy and lactation, while the diet of a second group of 16 pregnant females was treated with a chemical. For each cluster (litter consisting of the pups born to a female rat), the number n of pups alive at 4 days and the number y of pups that survived at 31 day lactation period were recorded. The data are given as a fraction in Table 4.
In Table 4, the numerator is the number of dead pups during 21 days lactation period, and the numerator is the number who survived past 4 days. The purpose of the experiment was to determine if the chemical treatment significantly affects the survival rate among the pups. That is, we need to test the null hypothesis . The data are presented in a 2 × 2 format in Table 5.
and , giving relative risk .
, , , , and , and the 95% confidence interval on AR is: 0.33 < AR < 0.45.
3. Split Cluster Design
Split-cluster experiments are being used by investigators in health sciences when naturally occurring aggregates of individuals with nested subgroups may be assigned to different treatments. Cited examples include split mouth trials, in which a subject’s mouth is divided into two segments that are randomly assigned to different treatment groups. In other situation, randomization to treatment conditions may be possible at the person level within the cluster. In this case, when the treatment conditions are available within each cluster, the design is referred to as a multisite or split cluster design (SCD). The major attractiveness of this design is that it removes a large portion of the inter-subject variation from the estimate of treatment effect; and hence has the potential to require a lesser number of subjects than a parallel arm design with the same power. When the response variable of interest is binary, statistical methods developed to evaluate the effect of intervention depends on non-parametric methods, as shown in  .
In this section we present the data layout for the SCD (see Table 6) and derive the large sample variance of the AR as a measure of effect size.
Under a similar set up to that we presented in the previous section and with appropriate change in notations the random variables and will have the same beta-binomial distributions, but they are no- longer independent.
Table 4. Weil’s data: Mortality due to exposure is a two arms clinical trial.
Table 5. Weil’s data collapsed in a 2 ´ 2 table.
Table 6. Data layout for the split-cluster design.
The correlation parameters are estimated as shown in Equations (7)-(9).
Although the AR estimator maintains the same expression under split clusters, its variance is affected by the correlations within the sub-clusters, and between units in the exposed and the non-exposed sub-clusters.
Using the delta method, we can therefore show that
where , , and
Here, is the intraclass correlation among the individuals in the sub-clus- ters of exposed, and is the intraclass correlation among the individuals in the sub-clusters of unexposed. Both correlations are estimated from the one-way ANOVA layout as explained in Equations (7)-(9). The cross-clusters correlation which is interpreted as an intercluster correlation denoted by is similarly estimated from the data by first ignoring the splitting structure of the data, and then use the one-way ANOVA to obtain the within and between mean squares. Substituting these quantities in (7) we obtain a moment estimator of .
Example 3: Split-Mouth Trial
For illustrating the proposed methodology, as a third example, we consider data from a split-mouth trial on 23 patients evaluating the effect of chlorhexidine in the treatment of gingivitis  . The data are presented in Table 7. The chlorhexidine and control treatments were randomly applied to four sites located in the patient’s left and right sides of the upper and lower jaws. We are interested here in testing the effect of treatment on the presence or absence of plaque, as based on the measurements taken two weeks after baseline and summarized in Table 8. The sample estimates and standard errors (SE) of P and Q, the proportion of patients having plaque in the chlorhexidine and control groups, are estimated at 0.89 (SE = 0.0343), and 0.77 (SE = 0.0491), respectively. The intra-class and inter-class correlation coefficients are estimated as , , as shown in the previous section, pooled estimate The sample estimate of the relative risk RR is 1.155.
, , .
, and the 95% CI on AR is (−0.112, 0.225).
4. Testing the Equality of Two Correlated AR Parameters
Interest is focused on studying the change in disease-exposure etiology under va- rying conditions. We illustrate this situation using the published data  .
For example in the case of family data we may be interested in evaluating the effect of disease status of a parental exposure variable on their siblings, which can be divided into males and females within the same family. In this case, we
Table 8. Disease distribution among males and females according to father (exposure variable) disease status.
have two correlated attributable risk estimator, one describing the disease-ex- pose etiology for males, and the other for females. The main interest here is to compare the AR of males to that of females from the same sib-ship.
Example 4: Correlated AR’s from Cross Sectional Study: Family Data
We now consider a highly structured clustered familial data that has a two level hierarchy with blood measurements taken on parents (level two) and their offspring (level one) together with other anthropometric features  . Familial data sets are known to have considerable “within-cluster” correlation due to the homogeneous nature of family members. The goal is to classify the offspring blood pressure status based on parents BP and other anthropometric features. The data set contains 223 families with a mean number of siblings equal to 3 siblings per family. The outcome variable in this data set is a binary variable defined as offspring blood pressure status. If simultaneously SBP > 130 and DBP > 80, then an offspring is considered diseased otherwise normal . The exposure variable in this example is whether a parent (here we select the father) has the condition (presence of exposure) or does not have the condition (absence of exposure). The data are presented in Table 8.
We present the general methodology as follows: Testing for gender difference in the population is formulated as testing the null hypothesis against a general unspecified alternative . Note that testing this null hypothesis is equivalent to testing .
Let the point estimators be denoted by and . The difference is asymptotically unbiased and has variance .
Hence the null hypothesis is rejected whenever falls in the interval or , where is the cut off point on the standard normal curve. With a slight difference in notation, is similar to the expression in Equation (6). We derive using the delta method. In general, the data will have a structure similar to that given in Table 9.
We define the moment estimator as before:
Let , then similar to the first situation, we have:
Table 9. Collapsed data for the analysis of correlated AR parameters.
, , , ,
and , , and
Moreover, is the intracluster correlation of the exposed clusters under condition, and is the intracluster correlation of the unexposed clusters under condition. They two parameters are estimated as described in (7).
For simplicity we assume that these correlations are constant among the exposed and unexposed.
Using the delta method we can show after some algebra that:
The correlation which, under both conditions is the average correlation among the responses, is estimated as described in Section 3.
The values inside the square bracket are given by:
Using the data in Table 8 we get:
Males: , , .
Females: , ,
, , and .
Therefore there is not enough evidence in the data to support the hypothesis of presence of gender differences for the paternal effect on the siblings’ hypertension status.
We carried out a Monte-Carlo study generating the observations from bivariate beta binomial distribution. We restricted our simulations to the situation when the intracluster and the cross clusters correlation are equal. We also assumed a fixed number of observations within each cluster. The purpose was to limit the number of scenarios under which we examine the properties of the proposed test statistic Z. The statistic is computed when both and are strictly positive with additional restriction, . If these conditions are not satisfied, the sample is replaced until a total of 1000 iterations are obtained for each parameter combination. Table 10 shows the empirical levels and powers of the test assuming that the number of clusters is the same under both conditions and for the exposed and the non-exposed clusters. The main conclusions from Table 10, is that the proposed test statistic hold its empirical Type I error rate levels. There is an increase in the test power when the correlation increases, and when the prevalence of exposure parameters and are away from the boundaries of the interval (0, 1). We also note an appreciate increase in the power when the number of clusters is above 50, and naturally when the tested parameters are well separated.
The population Attributable risk, like the odds ratio and relative risk is a measure of disease risk association. However it has a special appeal to public health epidemiologists as it measures the percent reduction in the chances of having the outcome among subjects who are exposed to the risk factor. Clearly, not everyone in the population is exposed to the risk factor. For example, in evaluating the relationship between consanguinity and the risk of PDA, not all parents are relatives. We assume say that 55% of women in the population (as in the Saudi traditional society) are married to a first cousin. To determine how much of a reduction there would be in PDA among CHD newborns we have 0.55 × 0.06 = 3.3%.
We have developed estimators of the variance and the confidence interval on AR when the units of sampling are aggregates of individuals under three study designs. In all situations the estimation of the intraclass correlation is crucial to
Table 10. Empirical type i error rates and powers based on 1000 replications from the bivariate beta binomial distribution. We set AR1 = 0.05, and therefore, .
conduct valid statistical inferences.
One of the objectives of this paper was to develop and evaluate simple test statistic that could be used to compare dependent attributable risks in the case of clustered dichotomous outcome variables.
1) Through simulations, a major finding of our work is that to test the equality of correlated attributable risks, either from cross sectional or cohort studies, one needs a much larger number of clusters than that expected to achieve high power.
2) An interesting extension of our study is to construct model-based inference on the AR. This would require the development of a semi parametric model similar to the generalized estimating equation, or a full probabilistic model such as generalized linear mixed model where the effect of multiple covariates may be accounted for.
3) A limitation of the simulation study is the restrictions that the number of observations in all clusters is held constant (balanced design) and that the within cluster and the cross cluster correlations are equal. The reason for this assumption is to limit the number of factors which affect the power so that reasonable conclusions can be made. But we believe that these restrictions should not affect the overall conclusions.
Conflict of Interest
The authors declare that there is no conflict of interest.