Estimation of MZ and DZ Twinning Rates

The estimation is usually based on WDR, assuming that the probabilities of a male birth and female birth are equal and that the sexes within a twin set are independent. Although the probability of a male birth is greater than .5, the reliability of the WDR's assumptions has never been conclusively verified or rejected (Fellman & Eriksson, Reference Fellman and Eriksson2006). According to WDR, the total number of DZ twin maternities is twice the number of twin maternities with opposite-sex (OS) twin sets. The number of MZ twin sets is the difference between the number of same-sex (SS) and OS twin sets.

Let the total number of maternities be n, the number of SS twin sets be n_S, and the number of OS twin sets be n_O. One can assume that the distribution of n_S, n_O, and n − n_s − n_o has a trinomial distribution. Let the probability for an OS twin maternity be p_O, for an SS twin maternity p_S and, consequently, the probability for other maternities 1 − p_S − p_O. Applying WDR, the MZ rate is m = p_S − p_O and the DZ rate is d = 2p_O. According to the maximum likelihood theory, the estimators of transformed parameters are the corresponding transformations of the estimators of the initial parameters. Assuming that WDR holds, we can estimate p_O and p_S and use the formulae $\hat m = \hat p_S - \hat p_O$ and $\hat d = 2\hat p_O$. This method does not solve if the proposed estimators $\hat m$ and $\hat d$ are unbiased estimators of the MZ and DZ rates, respectively. The unbiasedness depends on the reliability of Weinberg's assumptions. The likelihood function is

(1)\begin{equation} L(p_S ,p_O ) = \left( {1 - p_S - p_O } \right)^{n - n_S - n_O } p_S^{n_S } p_O^{n_O } .\end{equation}

The log-likelihood function is

\begin{eqnarray*} l(p_S ,p_O ) &=& (n - n_S - n_O )\ln (1 - p_S - p_O ) + n_S \ln (p_S )\nonumber\\ && +\, n_O \ln (p_O ). \end{eqnarray*}

After differentiation, we obtain

(2)\begin{eqnarray} \frac{{\partial l}}{{\partial p_S }} &=& - \frac{{n - n_S - n_O }}{{1 - p_S - p_O }} + \frac{{n_S }}{{p_S }}\,\quad {\rm and}\nonumber\\ \frac{{\partial l}}{{\partial p_O }} &=& - \frac{{n - n_S - n_O }}{{1 - p_S - p_O }} + \frac{{n_O }}{{p_O }}.\end{eqnarray}

The maximum of l(p_S, p_O) yields the estimators $\hat p_S = \frac{{n_S }}{n}$ and $\hat p_O = \frac{{n_O }}{n}$, being unbiased, efficient, and asymptotically normal. To obtain the covariance matrix, we differentiate the formulae in (2) once more and obtain

\begin{eqnarray*} \frac{{\partial ^2 l}}{{\partial p_S^2 }} &=& - \frac{{n - n_S - n_O }}{{\left( {1 - p_S - p_O } \right)^2 }} - \frac{{n_S }}{{p_S^2 }},\quad\nonumber\\ \frac{{\partial ^2 l}}{{\partial p_O^2 }} &=& - \frac{{n - n_S - n_O }}{{\left( {1 - p_S - p_O } \right)^2 }} - \frac{{n_O }}{{p_O^2 }}\,,\quad {\rm and}\quad \,\nonumber\\ \frac{{\partial ^2 l}}{{\partial p_O \partial p_S }} &=& - \frac{{n - n_S - n_O }}{{\left( {1 - p_S - p_O } \right)^2 }}. \end{eqnarray*}

Now, E(n_S) = np_S and E(n_O) = np_O and, consequently,

\begin{equation*} - E\left( {\frac{{\partial ^2 l}}{{\partial p_S^2 }}} \right) = \frac{{n\left( {1 - p_O } \right)}}{{p_S \left( {1 - p_S - p_O } \right)}}, \end{equation*}

\begin{equation*} - E\left( {\frac{{\partial ^2 l}}{{\partial p_O^2 }}} \right) = \frac{{n\left( {1 - p_S } \right)}}{{p_O \left( {1 - p_S - p_O } \right)}}\,, \end{equation*}

and

\begin{equation*} - E\left( {\frac{{\partial ^2 l}}{{\partial p_O \partial p_S }}} \right) = \frac{n}{{\left( {1 - p_S - p_O } \right)}}. \end{equation*}

The information matrix for the estimators $\hat p_S$ and $\hat p_O$ is

(3)\begin{eqnarray} \bm{\rm I}\left( {\hat{p}_{S},\hat{p}_{O} } \right) = \left[{\begin{array}{cc} {\displaystyle\frac{{n\left( {1 - p_{O} } \right)}}{{p_{S} \left( {1 - p_{S} - p_{O} } \right)}}} & {\displaystyle\frac{n}{{\left( {1 - p_{S} - p_{O} } \right)}}} \\[11pt] {\displaystyle\frac{n}{{\left( {1 - p_{S} - p_{O} } \right)}}} & {\displaystyle\frac{{n\left( {1 - p_{S} } \right)}}{{p_{O} \left( {1 - p_{S} - p_{O} } \right)}}} \\ \end{array}} \right]\nonumber\\ \end{eqnarray}

and the covariance matrix is

(4)\begin{eqnarray} &&\bm{\rm C}\left( {\hat{p}_S,\hat{p}_O } \right) = \bm{\rm I}^{ - 1} \left( {\hat{p}_S ,\hat{p}_O } \right)\nonumber\\ &&\quad = \frac{1}{n}\left[ {\begin{array}{cc} {\displaystyle p_S \left( {1 - p_{S} } \right)} & { - p_{S} p_{O} } \\[4pt] {\displaystyle - p_{S} p_{O} } & {p_{O} \left( {1 - p_{O} } \right)} \end{array}}\right].\end{eqnarray}

The estimators $\hat p_S$ and $\hat p_O$ have the variances ${\rm Var}(\hat p_S ) = \frac{{p_S (1 - p_S )}}{n}$ and ${\rm Var}(\hat p_O ) = \frac{{p_O (1 - p_O )}}{n}$ and the correlation coefficient ${\rm Cor}\left( {\hat p_S ,\hat p_O } \right) = - \sqrt {\frac{{p_S p_O }}{{\left( {1 - p_S } \right)\left( {1 - p_O } \right)}}}$.

The linear transformations $\hat m = \hat p_S - \hat p_O$ and $\hat d = 2\hat p_O$ are efficient and asymptotically normal and, if WDR holds, also unbiased. The estimators of these parameters are $\hat m = \hat p_S - \hat p_O = \frac{{n_S - n_O }}{n}$ and $\hat d = 2\hat p_O = \frac{{2n_O }}{n}$. The covariance matrix is

(5)\begin{eqnarray} \begin{array}{l} \bm{\rm C}(\hat m,\hat d)\\ \quad = \left[\! {\begin{array}{*{20}c} 1 & { - 1} \\ 0 & 2 \end{array}}\! \right]\left[\! {\begin{array}{*{20}c} \hfill {\displaystyle\frac{{p_S \left( {1 - p_S } \right)}}{n}} & \hfill {\displaystyle - \frac{{p_S p_O }}{n}} \\[6pt] \hfill {\displaystyle - \frac{{p_S p_O }}{n}} & \hfill {\displaystyle\frac{{p_O \left( {1 - p_O } \right)}}{n}} \end{array}} \!\right]\left[\! {\begin{array}{*{20}c} 1 & 0 \\ { - 1} & 2 \end{array}}\! \right] = \\[21pt] \quad = \left[\!\! {\begin{array}{*{20}c} \hfill {\displaystyle\frac{{p_S + p_O - \left( {p_S - p_O } \right)^2 }}{n}} & \hfill {\displaystyle - \frac{{2p_O \left( {1 + p_S - p_O } \right)}}{n}} \\[6pt] \hfill {\displaystyle - \frac{{2p_O \left( {1 + p_S - p_O } \right)}}{n}} & \hfill {\displaystyle\frac{{4p_O \left( {1 - p_O } \right)}}{n}} \end{array}}\!\! \right]. \\ \end{array}\nonumber\\[-8pt] \end{eqnarray}

We rewrite (5) by using the parameters m and d and obtain

(6)\begin{eqnarray} \bm{\rm C}(\hat{m},\hat{d}) = \frac{1}{n}\left[ {\begin{array}{cc} \hfill {m(1 - m) + d} & \hfill { - d(1 + m)} \\ \hfill { - d(1 + m)} & \hfill {d(1 - d) + d} \end{array}} \right].\end{eqnarray}

We see that

(7)\begin{eqnarray} {\rm Var}(\hat m) &=& \frac{{m(1 - m) + d}}{n},\quad\nonumber\\ {\rm Var}(\hat d) &=& \frac{{d(1 - d) + d}}{n}\,,\quad {\rm and}\quad\nonumber\\ {\rm Cov}(\hat m,\hat d) &=& - \frac{{d(1 + m)}}{n}.\end{eqnarray}

Bulmer (Reference Bulmer1970) presented slightly different approximate variance formulae. He assumed that the numbers of SS and OS twin sets are Poisson distributed but ignored that the numbers of OS and SS twin sets are correlated. His formulae are ${\rm Var}(\hat m) = \frac{{m + d}}{n}$ and ${\rm Var}(\hat d) = \frac{{2d}}{n}$, yielding figures that, compared with ours, are slightly too large.

If $\hat m$ and $\hat d$ are considered to be relative frequencies, one obtains the formulae ${\rm Var}(\hat m) = \frac{{m(1 - m)}}{n}$ and ${\rm Var}(\hat d) = \frac{{d(1 - d)}}{n}$, which are incorrect. The differences between the correct and incorrect formulae are positive and of the same order of magnitude (n ^{− 1}) as the correct ones and, consequently, the simple formulae should not be used, not even when the data sets are large. In addition, the erroneous variances are too small and may yield quite misleading test results and confidence intervals (CIs). For $\hat d$, the ratio between the correct variance estimate and the erroneous one is $\frac{{d(1 - d) + d}}{{d(1 - d)}} = \frac{{2 - d}}{{1 - d}} \approx 2$, and for $\hat m$ it is $\frac{{m(1 - m) + d}}{{m(1 - m)}} \approx \frac{{m + d}}{m}$. The latter ratio is approximately the ratio between the total and the MZ twinning rate and can be considerably larger than the ratio for $\hat d$.

Testing

According to (7), the estimators are negatively correlated and the correlation is usually rather strong. Therefore, the simultaneous testing of both MZ and DZ rates is to be preferred and a joint confidence region is better than two separate CIs. The information matrix is the inverse of the covariance matrix formula (6), which is

\begin{eqnarray*} &&\bm{\rm I}(\hat{m},\hat{d}) \nonumber\\ &&\quad = n\left[ {\begin{array}{cc} \hfill {\frac{{d(1 - d) + d}}{{d\left( {2m - 3md - 2m^{2} + d - d^2 } \right)}}} & \hfill {\frac{{d(1 + m)}}{{d\left( {2m - 3md - 2m^{2} + d - d^{2} } \right)}}} \\[8pt] \hfill {\frac{{d(1 + m)}}{{d\left( {2m - 3md - 2m^{2} + d - d^{2} } \right)}}} & \hfill {\frac{{m(1 - m) + d}}{{d\left( {2m - 3md - 2m^{2} + d - d^{2} } \right)}}} \end{array}} \right]\nonumber\\ &&\quad = \left[ {\begin{array}{cc} {I_{11} } & {I_{12} } \\ {I_{21} } & {I_{22} } \end{array}} \right]. \end{eqnarray*}

We construct the test statistic

(8)\begin{eqnarray} {\rm \chi }^{2} &=& \left[ {((\hat{m} - m)(\hat{d} - d)} \right] \bm{\rm I}\left( {\hat{m},\hat{d}} \right)\left[ {\begin{array}{c} {\hat{m} - m} \\ {\hat{d} - d} \end{array}} \right]\nonumber\\ &=& I_{11} \left( {\hat{m} - m} \right)^2 + 2I_{12} \left( {\hat{m} - m} \right)\left( {\hat{d} - d} \right)\nonumber\\ &&+\, I_{22} \left( {\hat{d} - d} \right)^{2} ,\end{eqnarray}

which is approximately χ² distributed with 2 degrees of freedom. If we want to test the given hypothetical rates m and d, we observe $\hat m$ and $\hat d$ and use the formula (8). It can also be used for the construction of a simultaneous confidence region for m and d. We observe that $P( {I_{11} ( {\hat m - m} )^2 + 2I_{12} ( {\hat m - m} )( {\hat d - d} ) + I_{22} ( {\hat d - d} )^2 \!\! \le\! {\rm \chi }_\alpha ^2 } )\break \ge 1 - {\rm \alpha }$ is a confidence region for the parameters (m, d) with the confidence level 1 − α. The obtained region is an ellipse in the (m, d) plane (cf. Figure 1).

FIGURE 1 Confidence regions for MZ and DZ rates for Åland–Åboland, 1653–1949. Using ${\rm SE}(\hat m)$ and ${\rm SE}(\hat d)$, we can construct the 95% CIs for the MZ rate and the DZ rate. For the Åland–Åboland data, we get the interval (2.16, 3.45) for m and (15.64, 17.32) for d, and the rectangle is constructed according to the individual CIs for m and d (Fellman & Eriksson, Reference Fellman and Eriksson2006).

For tests of a set of MZ and DZ rates, contingence tables cannot be used. An approximate test of a set of MZ or DZ rates can be performed in the following way. Consider MZ rates. Under the null hypothesis that the MZ rates are constant, being m for (t = 1, . . ., T), ${\rm Var}(\hat m_t ) = \frac{{m(1 - m) + d}}{{n_t }}$, where n_t is the sample size for t = 1, . . ., T. Define the weighted mean $\bar m = \sum_t {\frac{{n_t \hat m_t }}{{\sum {n_t } }}}$. Under the null hypothesis, $\bar m$ is the most efficient estimate of m. Consequently, ${\rm \chi }^2 = \frac{1}{c}\sum_{t = 1}^T {n_t \left( {m_t - \bar m} \right)^2 }$, where c = m(1 − m) + d, is asymptotically χ² distributed with T − 1 degrees of freedom. An approximate χ² test can be obtained if one uses $\hat c = \hat m(1 - \hat m) + \hat d$ as an estimate of c. For $\hat d$ one obtains a similar formula, but now c = d(1 − d) + d.

In numerical applications, the rates are usually given in per mille. This means that the theoretical variances should be scaled with a million and the estimators and standard errors with 1,000.