Some New Results on Overdispersed Models for the Human Sex Ratio Data

We have presented in this paper alternative computation to arrive at results obtained in Lindsey and Altham (1998) on the use of the beta-binomial, the double binomial and the multiplicative binomial models to the Human Sex Ratio data. Our procedure estimates fewer parameters and thus the goodness-of-fit test statistics computed here are based on more d.f. The results obtained, using SAS PROC NLMIXED agree with those presented in the Lindsey et al. paper. Further, we extend the study to two additional distributions, namely, the Com-Poisson binomial and the correlated binomial distributions. Our results give the estimated probabilities under all the models considered here and we are able to compute the success probability under the multiplicative binomial model for the global data.

Keywords: Beta-binomial; Multiplicative Binomial; Double Binomial; Com-Poisson; Correlated Binomial; Success Probability

Lindsey and Altham (1998) re-analyzed the human sex ratio data collected by Geissler (1889) [1,2]. The data relate to the distribution of the sexes of children in families in Saxony during 1876-1885. A reproduction of the data is also presented in Edwards (1958) for distribution of boys of size 2-13 [3]. We have in this paper further provided a presentation of the data arranged by family size and corresponding frequencies. The data has been arranged so that for instance the data used by Sokal and Rohlf (1969), p.80 can now be found in row F12 for n=12. Similarly, that used by Fisher (1958) for a family size of n=8 can be found in row F08 for n=8 [4]. The number of children in each family therefore varies from 1 to 12, and is denoted by n in this paper. A detailed description of the data can be found in Edwards (1958) [3].

Questions of interest from previous studies are succinctly stated in Lindsey and Altham and include ’whether the probability that the child is a boy varies among families and whether it can vary over time within a family’. As argued in Lindsey and Altham, because the Geissler study results in an aggregated data, both questions may not be distinguishable and that they will inherently manifest themselves as overdispersion with respect to the binomial model. As presented in Lindsey and Altham, we will look at whether variability changes with family size.

In this study, we shall compare the behaviors of the binomial model (BN) with those of five other models, namely, the beta-binomial (BB) Skellam (1948); the multiplicative binomial (MBM) Altham (1978), Elamir (2013); the Com-Poisson Binomial (CPM) Borges et al. (2013), the double binomial (DBM) Lindsey and Altham (1995) and the correlated binomial Kupper and Haseman (1978) distributions [5-10].

Lindsey & Altham (1998), employed the beta-binomial, the multiplicative and the double binomial models as overdispersed models for the data in Table 1. In their paper, the forms of the multiplicative and double binomial models employed were characterized with intractable normalizing constants c(n,ψ) and c(n,π) respectively. Consequently, these models were implemented in Lindsey and Altham (1998) by utilizing a generalized linear model with a Poisson distribution and log link to the frequency data. This approach has earlier being similarly employed in Lindsey and Altham (1995) [1,9]. This approach, which employs joint sufficient statistics in both distributions was proposed in Lindsey and Mersch (1992) [11]. Thus, both distributions are fitted using a Poisson regression model having sufficient statistics from both distributions as explanatory variables with the frequencies being the mean dependent variables. For the Double binomial model (DBM), the sufficient statistics are ylog(y) and (n-y)log(n-y). Similar sufficient statistics for the multiplicative binomial model (MBM) are y and y(n-y)with the offset being z=log(n y) for both models. However, this approach is not efficient. Their approach fits the following model:

resulting in,π_n 1/(1exp(-lr)). Here, we see that the logit of the probability of “success” is a linear function of family size n. This approach is very similar to those employed in Lefkopoulou et al. (1989), where,

The model in (1) is employed to model the probability of success for each of the distributions considered in this study, the exception being the Com-Poisson binomial model (CPM) described later in this paper [12].

Further, the dispersion parameters for each of the models are expressed as exp(a₀+na₁),n=2,3,....12. The exception here being the correlated binomial (CBM) and the Com-Poisson binomial, because the dispersion parameters θ_i and v_i in both models can take positive and negative values.

The formulation of the success probability model and the dispersion parameter models are as formulated in Lindsey and Altham (1998) [1]. However, because of the joint sufficient statistics’ approach employed in their paper, each model so formulated is based on thirteen parameters. In this study, the same models will be fitted, but with only four parameters, giving us more degrees of freedom for inference. Our models generate perfectly, the results in Lindsey and Altham (1998), even though we have employed only four parameters for each distribution [1]. We also examine the following variations of these models, for n=2,3,....12 family sizes.

(a) The model that assumes that for all n, we have constant probability of ’success’ and constant dispersion parameter (ϕ is that for the double binomial for instance). That is,

Since the data is binary, a binomial model applied to the entire data (excluding family size of one) gives an estimated probability of males ’success’ to be 0.5150. The dispersion parameter estimate from this model (X²/d.f.) is 9117.9513 which is very much greater than 1, indicating a very strong over-dispersion in the data. Thus, the binomial model would not fit these data and we thus consider in the following sections, the five over-dispersed probability distributions employed in this study.

Lovison (1998) proposed an alternative form of the two-parameter exponential family generalization of the binomial distribution first introduced by Altham (1978) which itself was based on the original Cox (1978) representation as [6,13,14]:

Where 0 < ψ < 1 and ω > 0. When ω =1 the distribution reduces to the binomial with π =ψ . If ω=1, ,n→∞ and 0ψ→, then nψ→μ and the MBD reduces to Poisson (μ).

Elamir (2013) presented an elegant characteristics of the multiplicative binomial distribution including its four central moments. His treatment includes generation of random data from the distribution as well as the likelihood profiles and several examples.

Following Elamir (2013), the probability π of success for the Bernoulli trial, can be computed from the following expression in (4) as [7]:

with π defined as in (4), ψ therefore can be defined as the probability of success weighted by the intra-units association measure ω which measures the dependence among the binary responses of the n units. Thus if ω=1, then π=ψ and we have independence among the units. However, if 1 ω≠, then, πψ≠and the units are not independent.

The mean and variance of the MBM are given respectively as:

The Beta-Binomial distribution Skellam (1948) has the probability distribution den-sity function [5]:

The probability density function for the Com-Poisson Binomial distribution is given by:

With π€ (0,1) and v€ R . If ν = 1, the model reduces to the binomial distribution and values of v>1 indicate underdispersion, while values of v< 1 similarly indicate overdispersion with respect to the binomial distribution

An approximation to the CPB distribution in the limit n → ∞ and with λ = n^ν p is given in Shmueli et al. (2005) [15]. Following Borges et al. (2013), if we let θ_n be defined as [8]:

The various properties of the CPB or the Com-Poisson have been presented in various papers Borges et al. (2013), Shmueli et al. (2005), and Kadane et al. (2006) applied the CPB to the number of killings in rural Norway [8,15,16]. In this study, with θ defined as θ_n=α₀+nα₁ for the family sizes, then the estimated probabilities are obtained as π_n=ˆθ_n/1+θ_n.

In Feirer et al. (2013), the double binomial distribution was presented, having the pdf form [17]:

For y = 0, 1, . . . , n, and following Feirer et al. (2013), its log-likelihood is presented in (16c).

A fetal death model proposed by Kupper and Haseman (1978), which is a generalization of the binomial model has the form [10]:

In (15), θ_i is the covariance between any two responses within the same litter. The model above allows for negative intra-litter correlation, as distinct from the other models considered that only allow for positive dispersion parameters.

Because of the bounds on θ_i for this model, it is modeled here as θ_i = b₀ + n_i b_i,i =2, 3, · · · ,1 2 rather than taking its exponent.

For a single observation, the log-likelihoods for the multiplicative binomial, beta- binomial, double binomial, the Com-Poisson binomial, and the correlated binomial are displayed in expressions (16a) to (16e) respectively.

Maximum-likelihood estimations of the above models are carried out with PROC NLMIXED in SAS, which minimizes the function -LL(y,Θ) over the parameter space numerically. The default integral approximations of the marginal likelihood in PROC NLMIXED is the Adaptive Gaussian Quadrature as defined in Pinheiro and Bates (1995) and the optimization algorithm adopted in this study is the Newton-Raphson iterative procedure (NEWRAP). Convergence is often a major problem here and the choice of starting values is very crucial [18]. For each of the cases considered here, we are able to achieve rapid convergence for all the distributions.

As discussed earlier, following Lindsey & Altham, we model the probability, πn as a function of family size n, viz logit(π_n)=β₀+nβ₁ Similarly, the corresponding dispersion parameters are also modeled as: log(ω)=α₀+nα₁. Thus, for the multiplicative binomial, for instance,ω=exp(lp) , and both are functions of family size n. Further, for each model, we are estimating only four parameters where for instance, for the beta-binomial model for example, we are estimating only {α₀,α₁,β₀,β₁} as compared to 13 parameters that were estimated in Lindsey and Altham paper.

In Table 2 are presented the parameter estimates and their standard errors (in parentheses) under the five models for cases I to IV [19-21]. For the full models (Case IV), we see that only four parameters are estimated for each model in contrast to Lindsey and Altham (1998) thirteen parameters [1].

The -2LL values for the beta-binomial in case IV were 2174100.8 but we have reported this simply as 4100.8. That is, we are reporting only the last four digits. In each case, the values of 217 had been truncated for brevity only. What is obvious from Table 2, is that, irrespective of the model being employed, the data is better modeled with a variable dispersion parameter, that is, cases II and IV. In other words, models that assume uniform dispersion parameter (within a given distribution here) are not satisfactory for the data. The -2LL and goodness-of-fit test statistics, the Pearson’s

and the likelihood-ratio test

all indicate that the models in cases II and IV are much better. Further, of the two parsimonious models in cases II and IV, model IV which is equivalent to Lindsay and Altham (1998) model is most parsimonious. For instance, for the beta-binomial model, the difference in G2 for the two models is 125.5495 − 107.7504 =17.7991 on 1 d.f., (G² rather than X² is used here because of its partitioning property over X2) which clearly indicates that β1 is important in the model, given that α0, α1 and β0 are already in the model. Consequently, we would focus our attention on models in case IV for all the distributions.

For these five models, the best three models in case IV that best fit our data are in the order (a) Multiplicative binomial, (b) The beta-binomial and (c) the Com-Poisson binomial distributions. However, we present in Table 3, the estimated probabilities for each family size under each of these models.

In Table 3 are presented the estimated probabilities from each of these models. We present the computation of the first probabilities for family size n=6 as an example. For the beta-binomial, this equal π_i=1/[1+exp(-0.0498-0.0018*6)]= 0.5151.

Similarly, under the multiplicative binomial model,ψ^₆=1/[1+exp(-0.05219-(0.000864×6))]=0.51434.

However, this is not the required probability. The required probability,π¹₆ can be computed using (7) as:

We present in Table 4, the calculated values necessary to compute π₁, π₂ and the variance under the multiplicative binomial model. The table also includes computed estimated variances as well as the dispersion parameter ω. We shall refer to this table again later. Clearly, the beta-binomial, the double binomial and the correlated give estimated global family size probabilities that are very similar. Of the three models here, the beta-binomial fits best and we see that the estimated probabilities under the BB agree with those presented in Lindsey and Altham (1998), even though we have employed only four parameters to fit the models [1]. The Com-Poisson binomial seems to be underestimating the global probabilities slightly, but seems to fit the model better than the DBM and the CBM.

Similarly, the estimated π¹_n obtained (first probability) under the multiplicative binomial model are very consistent with those from the BB. The questions posed in Linsey and Altham regarding the good fit of the multiplicative binomial model, namely:

(a) Does it (i.e. the multiplicative model) have some biological significance similar to the idea for the beta-binomial model of the Bernoulli probability varying across the population according to a beta distribution?

(b) Can more be said about the parameters than the empirical statement that the probability and dispersion both increase with family size?

Based on our ability to now be able to compute the central moments of the multiplicative Binomial model (Elamir (2013)), the ψ in our model specification in (3) is not necessarily equivalent to the desired πn, our ability to now be able to compute the πⁱ_n,i=1,2,3,4 and in particular, the case when i = 1 which corresponds to our Bernoulli or BB parameter π_n from expressions in (7) answers the question poised in (a) [7]. Thus, π¹ is considered to be the ’probability of success while ψ is the probability of success weighted by the intra-association measure ω given the dependence among the binary responses of the units.

As to the second question, the results of ˆω_n in Table 5, indicate that ω^_n<1,∀i and decreases as the family size increases. Thus, while ω^₂=0.99050, that for ω^₁₂=0.97227. These indicate decreasing variability as family size increases. On the other hand, for the beta-binomial, both probability and dispersion parameter increase with family size. For values of n small, estimates of ω are very close to 1, indicating closeness to being independent observations and can be appropriately modeled with the binomial distribution.

Also, in Table 5, the computed dispersion parameters for the double binomial, ϕ and the Com-Poisson binomial, ν, indicate decreasing dispersion parameters as the family size increases. For both the beta and correlated binomial, the dispersion parameters however, increase with family size.

We present in this section results of goodness-of-fit tests for all the models considered in this study. From the global model for each distribution, and a given estimated set of parameters, the probabilities of Y=0,1 , . . . ,n and corresponding expected values are computed from the log-likelihood equations. For example, for the multiplicative model, this is, for a given family size n, the LL1 is as presented earlier in (16a). Hence, the corresponding probabilities are computed for Y=0,1 , . . . ,n as P^_Y=P(Y=y)=exp(LL1) and the corresponding expected values are ˆm_y=Np^_y, where N for the family of size 12 for instance is 6115. We give in Table 6, the computed probabilities, the expected values, and contributions to the likelihood-ratio test statistic G₂ and Pearson’s X₂ for the multiplicative binomial for n=2,3,4and n=12. Here, G² and X² are defined as

where ˆm_i=Np^_i= are the expected values as defined above for each family size n .

From Table 6, our computations for the G2 and X2 exclude the case when n =1, that is, family size of 1. Thus, we have the computations over 88 observations only. Note that for each n, the cumulative sum of probabilities (designated here by Σp_y) sums to 1 for each family size n. We are re-assured that for a given n these probabilities sum to one. We have presented the results for n=2,3,4 and n=12 only. The sample sizes for n = 2, 3, 4, 12 are respectively 179,892; 148,903; 120,137; and 6,115. The expected values under each N are provided under the column labeled m^_y . This approach is employed for all the models considered in this study and in Table 7 are the resulting GOFs for all the five distributions under case IV.

Table 7 gives the GOF for the five models (in comparison to the binomial model) and results here demonstrate that while the multiplicative binomial model is best, it is however not much better than the beta-binomial model. However, since the multiplicative model is not much difficult to model than the beta-binomial, it would be better to actually employ it. Our results in other studies regarding the two models suggest that the MBD may be a better choice. Clearly, both the BBD and the MBD show considerable improvements in model fits than the binomial model.

Our results here agree with those in Lindsey and Altham (1998) that the ”multi- plicative model fitted somewhat better than did the beta-binomial” [1]. Our procedure here estimates fewer parameters than those employed in Lindsey and Altham (1998) and the optimization techniques open to us with SAS PROC NLMIXED, makes the use of this procedure much easier to use [1]. The implementation and use of more models in this study gives a broader perspective on models used for data having binary outcomes, especially as are being applied to the Human Sex Ratio Data. The SAS programs employed for this analysis are available from the author.