Variance of Observed Proportions

In population studies of the relative frequency of twin maternities, the accepted variance formula for the observed rates is

(1) $$\begin{equation} {\rm{Var}}(\hat{p}) = \frac{{p(1 - p)}}{n}, \end{equation}$$

where p is the theoretical probability of success.

Consider n repeated Bernoulli trials, that is, trials where the probability of success in a specific trial is independent of earlier outcomes. Let S be the total number of successes (e.g., twin maternities) and let $\hat{p} = \frac{S}{n}$ be the observed proportion of successes. For a large n, both S and $\hat{p}$ are asymptotically normal. If the assumption holds that (a) the repetitions are independent and (b) the probability of success, p, is constant during the repetitions, then the variance formula (1) for the observed proportion $\hat{p}$ holds.

Usually there is no discussion about the premises, and so we can use formula (1). In some situations, however, such a discussion may be necessary. Let us consider the outcome of the maternities of a certain group of mothers. The maternities (the mothers) are the repetitions, and success is the birth of a twin set. If we consider different mothers, then we can assume independent repetitions. However, the constancy of the binomial proportion is more difficult to accept. We know that the probability of a twin maternity varies greatly, depending on several factors, particularly the individual factors of maternal age, parity and ethnicity. Hence, the variance formula (1) for the total TWR cannot hold exactly. Does it hold approximately and with what accuracy? Cramér (Reference Cramér1951, p. 206) has studied the effect of defects in assumption (b) on the accuracy of the variance formula.

In order to generate classes that are as homogeneous as possible, assume that the mothers are grouped in classes according to presumptive influential factors. The statistical analysis performed is based on the crucial assumption that the grouping is chosen before the outcome is known. Let the number of mothers in group number r be n_r (r = 1, . . ., R) and the total number of mothers be n, that is, n = ∑ _r n_r . Furthermore, let S_r (r = 1, . . ., R) be the observed number, and ${\hat{p}_r} = \frac{{{S_r}}}{{{n_r}}}$ be the observed TWR within group number r and let S = ∑ _r S_r . Assuming constant probability within the groups, the theoretical probability of twin maternities in the total population is $p = \frac{1}{n}\sum_r {{n_r}{p_r}} $ , where p_r is the group-specific probability. The total observed relative frequency of successes is

(2) $$\begin{equation} \hat{p} = \frac{1}{n}\sum\limits_r {{n_r}{{\hat{p}}_r}} = \frac{1}{n}\sum\limits_r {{S_r}} = \frac{S}{n}. \end{equation}$$

The expectation is $E(\hat{p}) = \frac{1}{n}\sum_r {{n_r}E({{\hat{p}}_r}) = \frac{1}{n}\sum_r}$\break ${ {{n_r}{p_r} = p} } $ , and $\hat{p}$ is an unbiased estimator of p. Note that the estimator $\hat{p}$ in (2) is the same, irrespective of whether we consider grouping or not. If we assume a constant probability within the classes but not between them, we obtain (Cramér, Reference Cramér1951: Fellman & Eriksson, Reference Fellman and Eriksson2004)

(3) $$\begin{eqnarray} {\rm{Var}}(\hat{p}) &=& \frac{1}{{{n^2}}}\sum\limits_r {n_r^2{\rm{Var}}({{\hat{p}}_r})} = \frac{1}{{{n^2}}}\sum\limits_r {n_r^2\frac{{{p_r}(1 - {p_r})}}{{{n_r}}}}\nonumber\\ &=& \frac{{p(1 - p)}}{n} - \frac{1}{{{n^2}}}\sum\limits_r {{n_r}{{({p_r} - p)}^2}} . \end{eqnarray}$$

This fundamental result can also be written

(4) $$\begin{equation} \frac{{p(1 - p)}}{n} = \frac{1}{{{n^2}}}\sum\limits_r {{n_r}{p_r}(1 - {p_r})} + \frac{1}{{{n^2}}}\sum\limits_r {{n_r}{{({p_r} - p)}^2}} . \end{equation}$$

The left-hand-side in (4) is the total variation, giving the variance when one ignores any grouping. The first part on the right-hand side is the variation within the groups and, according to (3), it is ${\rm{Var}}(\hat{p})$ when the grouping is considered. The second part is the variation between the groups. If, and only if, the rate is the same within all the groups, that is, p_r = p for all r = 1, . . ., R, the second term is equal to zero and all the variation is intra-group variation. In this case, the grouping factors have no influence on the incidence of twinning and the standard formula (1) for the variance holds. The first term on the right-hand side of (4) is zero if, and only if, p_r = 1 or p_r = 0. This means that the set of mothers is divided into groups consisting entirely of mothers with or without twin maternities. Such a grouping in surely homogeneous classes prior to the outcomes is possible only if every mother forms her own group. In this case, the total variation consists of inter-group variation and formula (1) holds.

The results obtained have an interesting interpretation. Without any grouping, the variance formula (1) holds. A grouping giving additional information about the variation of the probability reduces the variance. The reduction from (1) to ${\rm{Var}}(\hat{p}) = \frac{1}{{{n^2}}}\sum_r {{n_r}{p_r}(1 - {p_r})} $ indicates the effect of grouping. If the variation in the p_rs is large, then the grouping factors (race, maternal age, parity, marital status, time, rural or urban population, season, etc.) are informative with respect to the TWR. On the other hand, if the variation in the p_rs is very small then the grouping factors are of small informative value, and consequently, the reduction from (1) to (3) can be ignored. Furthermore, the overestimation indicates that all statistical tests give excessively low test values, that is, the tests are conservative.

According to Fellman and Eriksson (Reference Fellman and Eriksson2004), both the classical variance formula and the variance formula (3) are over-estimates as long as the groups are still heterogeneous. This follows from the fact that the formulae are based on the assumption that the probability of a twin maternity is constant within every group. Within these groups, additional factors, not considered, may cause unidentifiable heterogeneity. Such variations may be caused by inter-individual differences in the probabilities. Fellman and Eriksson (Reference Fellman and Eriksson2004) observed that the corrections are minute, but in studies where one considers small series, great differences between groups, or large numbers of groups (cf. grouping according to individual mothers, discussed above), the corrections may be notable. The formulae derived above are the only ones available for correcting the variance of group heterogeneity in the probabilities for twin maternities. It is very common that registered twinning data are pooled, and later it is impossible to split the data into more homogeneous groups.

Confidence Interval for an Unknown Proportion

Now, we consider confidence intervals (CIs) for an unknown proportion p. The ML estimator of p is $\hat{p} = \frac{S}{n}$ , where S is the number of successes (in our study, the number of twin maternities) and n is the total number of observations (all maternities). It is a well-known fact that $\hat{p}$ has the mean p, is asymptotically normal, and, if we consider ungrouped data or data for a specific group, the variance ${\rm{Var}}(\hat{p}) = \frac{{p(1 - p)}}{n}$ is given in (1). The standard method, at least in applied studies, is that the variance formula is estimated so that p is replaced by $\hat{p}$ , resulting in the approximate test statistic

(5) $$\begin{equation} z = \frac{{\hat{p} - p}}{{\sqrt {\frac{{\hat{p}(1 - \hat{p})}}{n}} }}. \end{equation}$$

Based on this formula, the 100(1-α)% CI, the so-called Wald's CI, is

(6) $$\begin{equation} \left( {\hat{p} - {z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }}\sqrt {\frac{{\hat{p}(1 - \hat{p})}}{n}} ,\;\hat{p} + {z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }}\sqrt {\frac{{\hat{p}(1 - \hat{p})}}{n}} } \right), \end{equation}$$

where z _1/2α is the (1-½ α) quantile for the standardized normal distribution.

Brown et al. (Reference Brown, Cai and DasGupta2001; Reference Brown, Cai and DasGupta2002) gave a thorough presentation of the problems concerning CIs for unknown proportions. They showed that the standard method (6) in universal use is riddled with problems. The actual coverage probability for p can differ significantly from the nominal confidence level at realistic and even larger than realistic sample sizes. The error comes from two sources: discreteness and skewness (p ≠ 0.5) in the underlying distribution. For a two-sided interval, the rounding error due to discreteness is asymptotically dominant, being of the order n ^−½. The error due to skewness is of the order n ⁻¹, but still important for even moderately large n. For one-sided intervals, the error caused by skewness can be larger than the rounding error. Brown et al. (Reference Brown, Cai and DasGupta2002) applied Edgewood expansions and found that two-term expansions fitted very well. This finding is a consequence of the fact that the first two terms in the Edgewood expansions are of the orders n ^−½ and n ⁻¹, respectively, being of the same order as the errors presented above.

The standard CI is centered at $\hat{p} = \frac{S}{n}$ . Although $\hat{p}$ is an unbiased ML estimator of p, as a center of a CI, it causes a systematic negative bias in the coverage. Brown et al. (Reference Brown, Cai and DasGupta2001) stressed that even in cases where the textbooks indicate that the formula (6) is safe, the coverage probability may differ markedly from the expected. As alternatives, they proposed the Wilson (Reference Wilson1927) or the Agresti–Coull CI (Agresti & Coull, Reference Agresti and Coull1998). The alternative methods are briefly presented below.

The Wilson Interval

Using the notations introduced above, we consider the exact test statistic

(7) $$\begin{equation} z = \frac{{\hat{p} - p}}{{\sqrt {\frac{{p(1 - p)}}{n}} }}, \end{equation}$$

which, on theoretical grounds, can be considered more reliable than the approximate one in (5). From (7) we obtain, after some calculations, the 100(1-α)% Wilson confidence interval

(8) \begin{eqnarray} \left( {\frac{{n\hat{p} + \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }^2}}{{n + z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }^2}} - \frac{{{z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }}\sqrt n }}{{n + z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }^2}}\sqrt {\hat{p}(1 - \hat{p}) + \frac{{z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }^2}}{{4n}}} ,\;\;}\right.\nonumber\\ \left.{\frac{{n\hat{p} + \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }^2}}{{n + z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }^2}} + \frac{{{z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }}\sqrt n }}{{n + z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }^2}}\sqrt {\hat{p}(1 - \hat{p}) + \frac{{z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }^2}}{{4n}}} } \right). \end{eqnarray}

The Agresti–Coull Interval

If we start from (8), use its mid-point $\tilde{p} = \frac{{n\hat{p} + \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }^2}}{{n + z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }^2}}$ as a modified rate, and calculate the interval analogously to (6), we then obtain the Agresti–Coull confidence interval

(9) $$\begin{equation} \left( {\tilde{p} - {z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }}\sqrt {\frac{{\tilde{p}(1 - \tilde{p})}}{{n + z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }^2}}} ,\;\;\tilde{p} + {z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }}\sqrt {\frac{{\tilde{p}(1 - \tilde{p})}}{{n + z_{\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} \alpha }^2}}} } \right). \end{equation}$$

This is related to the ‘adjusted Wald’ CI, defined according to the following advice for a 95% CI: ‘Add two successes and two failures and then use the Wald formula’.

Furthermore, their conclusion was that the Agresti–Coull interval dominates the other intervals in coverage but is, on average, slightly longer. The Wilson (Reference Wilson1927) interval is comparable with the Agresti–Coull interval in both coverage and length. In addition, Brown et al. (Reference Brown, Cai and DasGupta2001) stressed that if one takes simplicity of presentation and ease of computation into account, the Agresti–Coull interval, although a bit too long, could be recommended for use. In our opinion, the Wilson CI is easy to interpret and is comparable with the Agresti–Coull CI in simplicity; we therefore recommend both as improved alternatives to the standard CI.

In studies of multiple maternities, the choice of method is important, for the conditions are conflicting. In general, one can expect a large number of maternities (n), but the proportion of multiple maternities (p) is small, and consequently, the distribution is skewed. In addition, in this study of the TWRs in small isolates the number of maternities can be moderate and sometimes even small, and consequently, the choice of the CI formula is of great importance.