Estimation of Log-Linear-Binomial Distribution with Applications

Abstract

Log-linear-binomial distribution was introduced for describing the behavior of the sum of dependent Bernoulli random variables. The distribution is a generalization of binomial distribution that allows construction of a broad class of distributions. In this paper, we consider the problem of estimating the two parameters of log-linearbinomial distribution by moment and maximum likelihood methods. The distribution is used to fit genetic data and to obtain the sampling distribution of the sign test under dependence among trials.

1. Introduction

During the last three decades, a growing amount of literature has been observed in generalizing the classical discrete distributions. The main idea was to apply the extended versions of modeling different kinds of dependent count or frequency structure in various fields; see the work of Johnson et al. in [1], and of Bowman and George in [2], and of George and Bowman in [3], Yu and Zelterman [4, 5]. Failure to take account of correlation in the data will cause less precision for binomial-based estimates; see, for example, Kolev and Paiva [6].

As a generalization for the binomial distribution, Lovison [7] has derived the distribution of the sum of dependent Bernoulli random variables as an alternative of Altham's multiplicative-binomial distribution [8] from Cox's log-linear representation [9] for the joint distribution of 𝑛 binary-dependent responses and it will be called log-linear-binomial distribution. This distribution is characterized by two parameters and provides wider range of distributions than are provided by the binomial distribution where the log-linear binomial distribution includes underdispersion, overdispersion models and the binomial distribution as a special case.

In this paper, the methods of moment and maximum likelihood are used to estimate the two parameters of Lovison's log-linear binomial distribution. The variance-covariance matrix for the estimated parameters is obtained. Log-linear binomial distribution is also used to fit genetic data, and to obtain the sampling distribution of the sign test under dependence among trials.

In Section 2, we define Lovison's log-linear binomial distribution. Estimation of the parameters of the distribution based on moments and maximum likelihood methods are derived in Section 3. Ungrouped data is studied in Section 4 as a special case. Two applications are given in Section 5.

2. Lovison's Log-Linear-Binomial Distribution

Consider the random vector 𝐙=[𝑍₁,…,𝑍_𝑛]^′, 𝑍_𝑖 being a binary response which measures whether some event of interest is present, “1”, or absent, “0” for a sample of 𝑛 units and 𝑌_𝑛=∑^𝑛_𝑖=1𝑍_𝑖 denotes the sample frequency of successes. To accommodate the possible dependence among 𝑍_𝑖, and under the assumption that the units are exchangeable Lovison [7] has obtained the distribution of 𝑌_𝑛 as 𝑃𝑌_𝑛=𝑦=^𝑛𝑦𝜓^𝑦1−𝜓^𝑛−𝑦𝜔^{𝑦(𝑛−𝑦)}∑^𝑛_𝑡=0(^𝑛_𝑡)𝜓^𝑡1−𝜓^𝑛−𝑡𝜔^{𝑡(𝑛−𝑡)},𝑦=0,1,…,𝑛, (2.1) where 0<𝜓<1 and 𝜔>0 are the parameters; for more details about this distribution; see the work of Lovison in [7]. This distribution provides a wider range of distributions than is provided by the binomial distribution, for example, Figures 1 and 2 show the distribution of 𝑌_𝑛 for 𝑛=10, 𝜓=0.2,0.5, and different values of 𝜔. For the values of 𝜔>1, the distribution is sharper in the middle than the binomial.

(a)

(b)

(c)

(d)

Figure 1 

The distribution of 𝑌𝑛 for different values of 𝜔, 𝜓=0.2, and 𝑛=10.

(a)

(b)

(c)

(d)

Figure 2 

The distribution of 𝑌𝑛 for different values of 𝜔, 𝜓=0.5, and 𝑛=10.

As can be seen from the figures for some values of 𝜔<1, the distribution can be 𝑈, bimodal and unimodal shapes. The expected value and the variance of 𝑌_𝑛 are given by𝜇=𝐸𝑌_𝑛𝜎=𝑛𝜋,₂=𝑉𝑌_𝑛=𝑛𝜋(1−𝜋)+𝑛(𝑛−1)𝜋₁−𝜋². (2.2) Note that 𝑉𝑌_𝑛=𝑉𝑌_𝑏+𝑛(𝑛−1)Cov𝑍_𝑘,𝑍_ℎ (2.3) for 𝑘≠ℎ and Cov stands for covariance. The variance of the binomial is 𝑉𝑌_𝑏=𝑛𝜋(1−𝜋), (2.4) and the covariance of 𝑍_𝑘 and 𝑍_ℎ is Cov𝑍_𝑘,𝑍_ℎ=𝜋₁−𝜋², (2.5) where 𝜅𝜋=𝜓^𝑛−1𝜓,𝜔𝜅_𝑛,𝜋𝜓,𝜔₁=𝜓²𝜅_𝑛−2𝜓,𝜔𝜅_𝑛.𝜓,𝜔 (2.6) Therefore, the variance of the log-linear binomial is equal to the variance of binomial when Cov(𝑍_𝑘,𝑍_ℎ)=0, greater than the variance of binomial when Cov(𝑍_𝑘,𝑍_ℎ)>0, and less than the variance of binomial when Cov(𝑍_𝑘,𝑍_ℎ)<0.

The expected value and the variance of 𝑌_𝑛 are nonlinear on both 𝜓 and 𝜔. For example, the nonlinearity in the variance of 𝑌_𝑛 is depicted in Figures 3 and 4.

(a)

(b)

(c)

(d)

Figure 3 

Variance of 𝑌𝑛 for various values of 𝜔 at each value of 𝜓 and 𝑛=25.

(a)

(b)

(c)

(d)

Figure 4 

Variance of 𝑌𝑛 for various values of 𝜔 at each value of 𝜓 and 𝑛=25.

3. Estimation of the Parameters

Let a random sample of 𝑅 sets of 𝑛 trials each be available, the number of given 𝑦 successes being 𝑓_𝑦(𝑦=0,1,2,…,𝑛), and 𝑅=∑^𝑛_𝑦=0𝑓_𝑦.

3.1. Method of Moments

We can use the first two sample moments to find moment estimates for 𝜓 and 𝜔 as follows. The first sample moment is 𝑚₁=^𝑛_𝑦=0𝑦𝑓_𝑦𝑅,𝑀₂=^𝑛_𝑦=0𝑓_𝑟𝑦−𝑚₁²𝑅. (3.1) Equating these sample moments to the corresponding population moments we obtain the estimates of 𝜓 and 𝜔 by solving the two equations𝑚₁𝜅=𝑛𝜓^𝑛−1𝜓,𝜔𝜅_𝑛,𝑀𝜓,𝜔₂=𝑛(𝑛−1)𝜓²𝜅_𝑛−2𝜓,𝜔𝜅_𝑛𝜅𝜓,𝜔+𝑛𝜓^𝑛−1𝜓,𝜔𝜅_𝑛𝜓,𝜔−𝑛²𝜓²𝜅_𝑛−1𝜓,𝜔𝜅_𝑛𝜓,𝜔². (3.2) The solution of these two equations needs numerical methods. The numerical solution of these two equations may be found by using nonlinear equation solver (nleqslv) in statistical 𝑅-software. The author has a program written in 𝑅 for finding the moment estimates upon request.

3.2. Method of Maximum Likelihood

The method of maximum likelihood provides estimators that have a reasonable intuitive basis and many desirable statistical properties. The likelihood of the sample can be written as 𝐿𝜓,𝜔∣𝑛,𝑦,𝑓=𝑅!^𝑛_𝑦=0𝑃(𝑦;𝜓,𝜔)^𝑓𝑦𝑓_𝑦!. (3.3) Take the logarithm of this function𝑙=log𝑅!+^𝑛_𝑦=0𝑓_𝑦log𝑃𝑦;𝜓,𝜔−^𝑛_𝑦=0log𝑓_𝑦! (3.4) The first partial derivatives with respect to 𝜓 and 𝜔 are 𝜕𝑙=𝜕𝜓^𝑛_𝑦=0𝑓_𝑦𝜕/𝜕𝜓𝑃𝑦;𝜓,𝜔𝑃,𝑦;𝜓,𝜔𝜕𝑙=𝜕𝜔^𝑛_𝑦=0𝑓_𝑦(𝜕/𝜕𝜔)𝑃𝑦;𝜓,𝜔𝑃.𝑦;𝜓,𝜔 (3.5) Let 𝑞₀=𝑛𝑦𝜓^𝑦1−𝜓^{(𝑛−𝑦)}𝜔^{𝑦(𝑛−𝑦)},𝑞₁=^𝑛_𝑡=0𝑛𝑡𝜓^𝑡1−𝜓^{(𝑛−𝑡)}𝜔^{𝑡(𝑛−𝑡)}. (3.6) Hence,𝑃=𝑞𝑦;𝜓,𝜔⁰𝑞₁. (3.7)

3.2.1. Estimation of 𝜓 and 𝜔

After simplification the first partial derivative for 𝜓 and 𝜔 can be written as 𝜕𝑙=∑𝜕𝜓^𝑛𝑦=0𝑦𝑓_𝑦𝜓−∑^𝑛_𝑦=0𝑛−𝑦𝑓_𝑦−1−𝜓𝑅𝑞^′_1𝜓𝑞₁,𝜕𝑙=𝜕𝜔^𝑛_𝑦=0𝑦𝑛−𝑦𝑓_𝑦𝜔−𝑞^′_1𝜔𝑅𝑞₁. (3.8) where 𝑞^′_1𝜓=^𝑛_𝑡=0𝑛𝑡𝜓^𝑡1−𝜓^{(𝑛−𝑡)}𝜔^{𝑡(𝑛−𝑡)}𝑡𝜓−𝑛−𝑡,𝑞1−𝜓_′1𝜔=^𝑛_𝑡=0𝑛𝑡𝑡(𝑛−𝑡)𝜓^𝑡1−𝜓^{(𝑛−𝑡)}𝜔^{𝑡(𝑛−𝑡)−1}. (3.9) By solving 𝜕𝑙/𝜕𝜓=0 and 𝜕𝑙/𝜕𝜔=0 simultaneously we obtain the maximum likelihood estimates 𝜓_𝑚𝑙 and 𝜔_𝑚𝑙 under the sufficient condition 𝜕²𝑙/𝜕𝜓²<0 and 𝜕²𝑙/𝜕𝜔²<0.

The numerical solution for these two equations may be obtained using nonlinear equation solver (nleqslv) in 𝑅-software. The author has a program written in 𝑅 for finding the maximum likelihood estimates upon request.

3.2.2. The Asymptotic Variances

Suppose that 𝜃=[𝜓𝜔]^′ and under certain regularity conditions the information matrix is𝐼_𝑖𝑗(𝜃)𝜕=−𝐸²𝜕𝜃_𝑖𝜕𝜃_𝑗log𝐿(𝜃). (3.10) The asymptotic variance-covariance matrix can be obtained as 𝑉𝜃≈𝐼_𝑖𝑗(𝜃)⁻¹. (3.11) To find the information matrix for the log-linear binomial distribution, we note that 𝜕²𝑙𝜕𝜓²∑=−^𝑛𝑦=0𝑦𝑓_𝑦𝜓²−∑^𝑛_𝑦=0𝑛−𝑦𝑓_𝑦1−𝜓²−𝑞₁𝑞^′′_1𝜓𝑅−𝑞^′2_1𝜓𝑅𝑞²₁,𝜕²𝑙𝜕𝜔²=−^𝑛_𝑦=0𝑦(𝑛−𝑦)𝑓_𝑦𝜔²−𝑞₁𝑞^′′_1𝜔𝑅−𝑞^′2_1𝜔𝑅𝑞²₁,𝜕𝜕𝜓𝜕𝑙𝜕𝜔=−𝑅𝑞₁𝑞^′′_1𝜔.𝜓−𝑅𝑞^′_1𝜔𝑞^′_1𝜓𝑞²₁, (3.12) where 𝑞^′′_1𝜓=^𝑛_𝑡=0𝑛𝑡𝜓^𝑡1−𝜓^{(𝑛−𝑡)}𝜔^{𝑡(𝑛−𝑡)}𝑡(𝑡−1)𝜓²−2𝑡(𝑛−𝑡)𝜓+()1−𝜓𝑛−𝑡)(𝑛−𝑡−11−𝜓²,𝑞^′′_1𝜔=^𝑛_𝑡=0𝑛𝑡𝑡(𝑛−𝑡)[𝑡(𝑛−𝑡)−1]𝜓^𝑡1−𝜓^{(𝑛−𝑡)}𝜔^{𝑡(𝑛−𝑡)−2},𝑞^′′_1𝜔⋅𝜓=^𝑛_𝑡=0𝑛𝑡𝑡(𝑛−𝑡)𝜓^𝑡1−𝜓^{(𝑛−𝑡)}𝜔^{𝑡(𝑛−𝑡)−1}𝑡𝜓−𝑛−𝑡.1−𝜓 (3.13) Taking the expectation, the information matrix is obtained as 𝐼(𝜃)=𝑅𝜇𝜓²+𝑅𝑛−𝜇1−𝜓²+𝑞₁𝑞^′′_1𝜓𝑅−𝑞^′2_1𝜓𝑅𝑞²₁𝑅𝑞₁𝑞^′′_1𝜔⋅𝜓−𝑅𝑞^′_1𝜔𝑞^′_1𝜓𝑞²₁𝑅𝑞₁𝑞^′′_1𝜔⋅𝜓−𝑅𝑞^′_1𝜔𝑞^′_1𝜓𝑞²₁𝑅𝜇𝑛−𝜇−𝜎²𝜔²+𝑞₁𝑞^′′_1𝜔𝑅−𝑞^′2_1𝜔𝑅𝑞²₁, (3.14) where 𝐸[∑^𝑛_𝑦=0𝑦(𝑛−𝑦)𝑓_𝑦]=𝑛𝐸(∑^𝑛_𝑦=0𝑦𝑓_𝑦)−𝐸(∑^𝑛_𝑦=0𝑦²𝑓_𝑦)=𝑛𝑅𝜇−𝑅(𝜇²+𝜎²). The variance-covariance matrix will be 𝑉𝜃≈[𝐼(𝜃)]⁻¹. (3.15) The estimated variance-covariance matrix is 𝑉𝜃≈𝐼(𝜃)⁻¹. (3.16) The author has a program written in 𝑅 for finding the estimated variance-covariance matrix upon request.

4. Special Case: Ungrouped Data

If the values of Bernoulli random variables 𝑍_𝑖 are known, the parameters 𝜓 and 𝜔 can be estimated as follows. By noticing that in a vector of binary responses 𝑧 there are 𝑛(𝑛−1)/2 pairs of responses, and if the order is irrelevant three types of pairs are distinguishable: (𝑛−𝑦)(𝑛−𝑦−1)/2 pairs of (𝑧_𝑘=0,𝑧_ℎ=0), 𝑦(𝑦−1)/2 pairs of (𝑧_𝑘=1,𝑧_ℎ=1), and (𝑛−𝑦)𝑦 pairs of (𝑧_𝑘=0,𝑧_ℎ=1), or (𝑧_𝑘=1,𝑧_ℎ=0)where 𝑦=∑^𝑛_𝑘=1𝑧_𝑘; see the work of Lovison in [7]. Therefore, the estimate of 𝜔 is 1𝜔=CPR. (4.1) The estimated cross-product ratio (CPR) is CPR=0.25𝑦𝑦−1𝑛−𝑦𝑛−𝑦−1#𝑧_𝑘=0,𝑧_ℎ=1#𝑧_𝑘=1,𝑧_ℎ=0. (4.2) To obtain the maximum likelihood estimate of 𝜓, we need to solve 𝑑𝑙=𝑦𝑑𝜓𝜓−𝑛−𝑦−𝑞1−𝜓𝑞^′1𝜓₁=0. (4.3) The estimated variance of 𝜓 can be obtained when 𝑅=1.

5. Applications

5.1. Sampling Distribution of the Sign Test for Comparing Paired Sample

The sign test is a nonparametric test which makes very few assumptions about the nature of the distributions under test. It is for use with two repeated (or correlated) measures, and measurement is assumed to be at least ordinal. The usual null hypothesis for this test is that there is no difference between the two treatments (groups, 𝑔_𝑘, 𝑘=1,2). Formally, let 𝜏=𝑃(𝑋>𝑌), and then test the null hypothesis  𝐻₀: 𝜏=0.5 for no differences against 𝐻₁: 𝜋≠0.5 for differences. The sign test can be written as 𝑠_𝑑=^𝑛_𝑖=1𝐼𝑔_1𝑖−𝑔_2𝑖=^𝑛_𝑖=1𝑍_𝑖=𝑌_𝑛, (5.1) where 𝑍_𝑖=1,𝑔_1𝑖−𝑔_2𝑖>0(+),0,𝑥𝑔_1𝑖−𝑔_2𝑖<0(−). (5.2) Under the assumptions of two outcomes, fixed probability of success, and independent trials it is assumed that the sampling distribution of 𝑠_𝑑 is binomial distribution. The rejection region is 𝑃_value=𝑃reject𝐻₀∣𝐻₀istrue,𝑛,𝑦_𝑑=^𝑦𝑑_𝑗=0𝑛𝑗0.5^𝑛 (5.3) and 𝑦_𝑑=min(#1,#0). For a two-tailed test we reject 𝐻₀ if 2𝑃_value≤𝛼, else we do not reject 𝐻₀, where 𝛼 is prespecified value.

5.1.1. Sign Test Under Dependence of the Trials

Suppose the assumption of mutual independence in the data is violated and the trails are dependent; see, for example, the work of Tallis in [10] and Luceño [11]. In this case, we suggest the log-linear binomial distribution as a sampling distribution of 𝑠_𝑑=∑^𝑛_𝑖=1𝑍_𝑖 rather than the binomial distribution. Let 𝜏=𝑃(𝑋>𝑌)=𝑃(𝑍=1)=𝜓(𝜅_𝑛−1(𝜓,𝜔)/𝜅_𝑛(𝜓,𝜔)) represent the probability of success. Then the null hypothesis 𝐻₀: 𝜏=0.5 is equivalent to 𝐻₀: 𝜓=0.5. Therefore, under the null hypothesis, the rejection region can be obtained as 𝑃_value⋅lb=𝑃reject𝐻₀∣𝐻₀istrue,𝑛,𝑦_𝑑,𝜔=^𝑦𝑑_𝑗=0^𝑛𝑗0.5^𝑛𝜔^{𝑗(𝑛−𝑗)}∑^𝑛_𝑡=0(^𝑛_𝑡)0.5^𝑛𝜔^{𝑡(𝑛−𝑡)}, (5.4) where 𝑦_𝑑=min(#1,#0) and 𝜔=1/√CPR. For a two-tailed test, we reject 𝐻₀ if 2𝑃_value.lb≤𝛼, else do not reject 𝐻₀ and 𝛼 is a prespecified value (𝐻₀: there is no difference against 𝐻₁: there is difference.)

Example 5.1. A physiologist wants to know if monkeys prefer stimulation of brain area A to stimulation of brain area B. In the experiment, 14 rhesus monkeys from the same family are taught to press two bars. When a light comes on, presses on Bar 1 result in stimulation of area A and presses on Bar 2 result in stimulation of area B. After learning to press the bars, the monkeys are tested for 15 minutes, during which time the frequencies for the two bars are recorded. The data are shown in Table 1.
Using the binomial distribution, we have 𝑛=14, 𝑦_𝑑=3, and 𝑃_value is2𝑃_value=³_𝑗=0𝑗140.5¹⁴=2×0.0287=0.0574. (5.5) Since 0.0574>0.05, we do not reject 𝐻₀.
Using the log-linear binomial distribution, we have 𝑛=14, 𝑦_𝑑=3, and1𝜔=(55)(3)/(25)(8)=1.1, (5.6) and the 𝑃_value is 2𝑃_value.lb=³_𝑗=0^𝑗140.5¹⁴1.1^{𝑗(14−𝑗)}∑¹⁴_𝑡=0𝑡¹⁴0.5¹⁴1.1^{𝑡(14−𝑡)}=2(0.00687)=0.0137, (5.7) where 0.0137<0.05, so we reject 𝐻₀. That is, we would conclude that monkeys prefer stimulation in brain area B to stimulation in area A. Note that the rejection of 𝐻₀ agrees with Wilcoxon Signed-Ranks test for the same data; see the work of Weaver in [12] and Siegel and Castellan in [13].

Table 1 

Number of bar presses in brain stimulation experiment.

Power of the Test
Following the method given by Groebner et al. in [14], we may use the normal approximation to study the power of the sign test at 𝑛=14, 𝜔=1.1 and using one-tail test for simplicity as follows. The power of the test is power=1−𝛽,𝛽=𝑃accept𝐻₀∣𝐻₀isfalse. (5.8) For one-tailed test 𝐻₀: 𝜏≤0.5 and 𝐻₁: 𝜏>0.5 using binomial distribution with 𝑛=14, 𝜏=0.5 and 𝜎=14(0.5)(0.5)=1.8708. The critical value with 𝛼=0.05 is 𝑦_0.05,𝑏=14(0.5)+1.645(1.8708)=10.077. (5.9) The power of the test is power_𝑏≈1−𝑃𝑍<10.077−14𝜏₁1.8708, (5.10) for 𝜏₁>0.5.
Using log-linear binomial distribution with 𝑛=14, 𝜓=0.5, and 𝜔=1.1, we obtain 𝜎=1.4616, then under 𝐻₀ and 𝛼=0.05, the critical value is𝑦_0.05,lb=14(0.5)+1.645(1.4616)=9.404. (5.11) The power of the test under log-linear binomial distribution is power_𝑙𝑏≈1−𝑃𝑍<9.404−14𝜏₁1.4616, (5.12)𝜏₁>0.5.
Power of the sign test is given in Table 2. In this case, the log-linear binomial distribution shows improvement in the power of the sign test over the sign test under binomial; for example, when 𝜏₁=0.70, the power increases from 0.44 to 0.61 about 1.40 times.

Table 2 

Power of the sign test under binomial and log-linear binomial distribution for one-tail test with 𝑛=14, 𝜔=1.1, and 𝛼=0.05 based on normal approximation.

5.2. Fitting Genetic Data

The data are taken from Salmaniyaa hospital records in Bahrain for a genetic study on the gender ratio. Table 3 shows the number of male children in 3475 families with 7 children.

Table 3 

Numbers of male children in 3475 families of size 7

The first two sample moments are 𝑚₁=3.12057,𝑀₂=1.21567. (5.13) If we use binomial distribution, the estimated value of 𝜏 from the sample is 𝜏=3.120577=0.4458. (5.14) This value of 𝜏 was used to obtain the expected frequencies shown in Table 3. The value of 𝜒²=193.01 with 7 degrees of freedom gives𝑝=0<0.05, the simple binomial model has to be rejected. If we use log-linear binomial distribution, and fit the data by maximum likelihood, we find that 𝜓=0.423,𝜔=1.14543. (5.15) The estimated variance-covariance matrix is 𝑉𝜃≈70974.89107162510716.258927.83⁻¹=0.000017−.000021−.0000210.000137. (5.16) The expected frequencies based on these estimates are shown in Table 3. The value of 𝜒²=1.288 with 6 degrees of freedom gives𝑝=0.97>0.05. Thus, the log-linear binomial distribution provides a good fit to the observed data.

6. Conclusion

The parameters of the log-linear binomial distribution were estimated by the moment and maximum likelihood methods. Both methods needed solving nonlinear equations to obtain the estimators of the parameters. We used nonlinear equation solver (nleqslv) package in statistical 𝑅-software to find the estimates of the parameters. The variance-covariance matrix for the maximum likelihood estimates was obtained. Moreover, the sampling distribution of the sign test was studied when trials are dependent. A set of genetic data from Salmaniyaa hospital in Bahrain has been fitted using log-linear binomial distribution. The fit is found preferable over fitting the binomial distribution.

Acknowledgment

The author wishes to thank the two referees for their valuable comments that enhanced the contents and presentation of the paper.