Note on information bias and efficiency of composite
likelihood

arXiv:1612.06967v1 [math.ST] 21 Dec 2016

Ximing Xu, Nancy Reid and Libai Xu

Abstract
Does the asymptotic variance of the maximum composite likelihood estimator of
a parameter of interest always decrease when the nuisance parameters are known?
Will a composite likelihood necessarily become more efficient by incorporating additional independent component likelihoods, or by using component likelihoods with
higher dimension? In this note we show through illustrative examples that the answer to both questions is no, and indeed the opposite direction might be observed.
The role of information bias is highlighted to understand the occurrence of these
paradoxical phenomenon.
Key words: Pairwise likelihood; estimating function; Bartlett’s second identity;
Godambe information matrix; nuisance parameter.

1

Introduction

The likelihood function for a complex multivariate model may not be available or very difficult to evaluate, and a composite likelihood function constructed from low-dimensional
marginal or conditional distributions has become a popular alternative (Varin, 2008;
Varin, Reid & Firth, 2011). Suppose Y is a p-dimensional random vector with prob1

ability density function f (y; θ), with a q-dimensional parameter vector θ ∈ Θ. Given
a set of likelihood functions Lk (θ; y), k = 1, . . . , K, defined by the joint or conditional
densities of some sub-vectors of Y, the composite likelihood function (Lindsay, 1988) is
defined as
CL(θ; y) =

K
Y

Lk (θ; y)wk ,

k=1

where the wk ’s are nonnegative weights and the component likelihood Lk (θ; y) might depend only on a sub-vector of θ. The choice of Lk (θ; y) and the weights {wk } is critical
for improving the efficiency of the resulting statistical inference (Lindsay, 1988; Joe &
Lee, 2009; Lindsay, Yi & Sun, 2011). In this paper we focus on the two most commonly used composite likelihood functions in literature, independence likelihood and
pairwise likelihood, which are defined as CLind (θ; y) =
Qp−1 Qp
r=1

s=r+1

Qp

r=1

f (yr ; θ ) and CLpair (θ; y) =

f (yr , ys ; θ ), respectively. Given a random sample {y(1) , . . . , y(n) }, where

each y(i) is a p-dimensional vector, the composite log-likelihood function is

c`(θ; y) =

n
X

(i)

c`(θ; y ) =

i=1

n
X

log CL(θ; y(i) ),

i=1

and the maximum composite likelihood estimator (MCLE) is θ̂CL = arg maxθ c`(θ; y).
In addition to the computational simplicity, the composite likelihood function has
many appealing theoretical properties. In particular, under some regularity conditions,
θ̂CL is consistent and asymptotically normally distributed with variance equal to the
inverse of the Godambe information matrix: G(θ) = H(θ)J −1 (θ)H(θ) (Lindsay, 1988;
Varin, 2008; Xu & Reid, 2011). Here H(θ) = E{−∇θ uc (θ; y)} is the sensitivity matrix,
and J(θ) = Varθ {uc (θ; y)} is the variability matrix, with the composite score function
2

uc (θ; y) = ∇θ c`(θ; y). Throughout this paper we use I(θ) to denote the Fisher information
matrix of the full likelihood function. Given two composite likelihood functions CL1 (θ; y)
and CL2 (θ; y), CL2 (θ; y) is said to be more efficient than CL1 (θ; y) if CL2 (θ; y) has a
greater Godambe information matrix than CL1 (θ; y) in the sense of matrix inequality. It
is well known that the full likelihood function is more efficient than any other composite likelihood function under regularity conditions (Godambe, 1960; Lindsay, 1988), i.e.
I(θ) − G(θ) is non-negative definite.
In general, the second Bartlett identity does not hold for composite likelihood functions, i.e. H(θ) 6= J(θ). After Lindsay (1982), we call a composite likelihood CL(θ; y)
information-unbiased if H(θ) = J(θ), and information-biased, otherwise. Composite
likelihood-based inferential tools have been developed for hypothesis testing (Chandler
& Bate, 2007; Pace, Salvan, & Sartori, 2011) and model selection (Varin & Vidoni, 2005;
Gao & Song, 2010). Information bias of a composite likelihood can make the resulting
inference more difficult. For example, if the composite likelihood is information-unbiased,
the likelihood ratio statistic has the same asymptotic chi-square distribution as its full
likelihood counterpart. On the other hand, if it is information-biased the likelihood ratio
statistic converges in distribution to a weighted sum of some independent χ2 random variables (Kent, 1982). Adjustments have been proposed to the information-biased composite
likelihood ratio statistic such that the adjusted statistic has an asymptotic chi-square distribution (e.g., Chandler & Bate, 2007; Pace, Salvan, & Sartori, 2011).
The full likelihood function is information-unbiased, but an information-unbiased composite likelihood is not necessarily fully efficient. In fact, any component likelihood function Lk (θ; y) is information-unbiased. More generally, any composite likelihood function
3

as the product of component likelihoods with mutually uncorrelated score functions is
information-unbiased. As an example, consider a p-dimensional vector Y = (Y1 , . . . , Yp )T
with density function f (y1 , . . . , yp ; θ), and defining f (y1 | y0 ; θ) = f (y1 ; θ). It is easy to
show that the covariance between the score function of f (yi | y1 , . . . , yi−1 ; θ) and the score
function of f (yj | y1 , . . . , yj−1 ; θ) is zero for any i 6= j ∈ {1, . . . , p}. Hence any composite
likelihood of the form
Y

f (yi | y1 , . . . , yi−1 ; θ),

i∈A

where A ⊆ {1, . . . , p}, is information-unbiased. Conversely, an information-biased composite likelihood function can be fully efficient. The pairwise likelihood function for the
equal-correlated multivariate normal model in Section 2 is fully efficient when estimating
the common variance σ 2 and the correlation coefficient ρ (Cox and Reid, 2004), but it
is not information-unbiased (Pace, Salvan, & Sartori, 2011). A sufficient and necessary
condition for a composite likelihood to be fully efficient is given in the following theorem.
Theorem 1.1 Theorem 1.Suppose the full likelihood function L(θ | y) has the score
function u(θ) and Fisher information I(θ). Then, for any composite likelihood function
CL(θ | y) with the score function uc (θ), sensitivity matrix H(θ), variability matrix J(θ)
and Godambe information G(θ), G(θ) = I(θ) if and only if u(θ) = H(θ)J(θ)−1 uc (θ) + b(θ)
with probability 1 for a constant vector b(θ) with respect to the random vector y.
Proof

ProofIt is easy to show that H(θ) =Cov(uc (θ), u(θ)) (Lindsay, 1988). As I(θ) =Var(u(θ))

and J(θ) =Var(uc (θ)), the result follows as the difference I(θ)−G(θ) = I(θ)−H(θ)J(θ)−1 H(θ)
is the covariance matrix of u(θ) − H(θ)J(θ)−1 uc (θ).
Theorem 1 with b(θ) = 0 gives a sufficient condition for the maximum composite
4

likelihood estimator to coincide with the MLE (Kenne Pagui, Salvan & Sartori, 2014),
which is satisfied by the pairwise likelihood for closed exponential family models (Mardia
et al., 2009). In particular, the equicorrelated multivariate normal model with unknown
variance σ 2 and correlation coefficient ρ belongs to the closed exponential family, and
it has been shown that b(θ) = 0 and u(θ) = H(θ)J(θ)−1 uc (θ) or equivalently, uc (θ) =
J(θ)H(θ)−1 u(θ) for the pairwise likelihood function (Pace, Salvan, & Sartori, 2011). Its
application in more general exponential family models has been studied in Kenne Pagui,
Salvan, & Sartori (2014b).
In this paper we explore the impact of information bias on the composite likelihood
inference in more detail. In Section 2 we show through the equicorrelated multivariate
normal model that an information-biased composite likelihood may lead to less efficient
estimates of the parameters of interest when the nuisance parameters are known. A sufficient condition is also provided for the occurrence of such a paradoxical phenomenon. We
would expect that a more efficient composite likelihood can be obtained by incorporating additional independent component likelihoods or using higher dimensional component
likelihoods. However such strategies do not always work for information-biased composite
likelihood functions, as shown in Section 3. We conclude with a discussion in Section 4.

2

Composite likelihood with known nuisance parameters

In the presence of nuisance parameters, it is well known that the maximum likelihood
estimator of the parameter of interest will have a smaller asymptotic variance when the
5

nuisance parameters are known. It is easy to check that this also holds for informationunbiased composite likelihood functions. Suppose the q-dimensional parameter vector θ
is partitioned as θ = (ψ, λ), where ψ is a q1 -dimensional parameter vector of interest and
λ is a q2 -dimensional nuisance parameter vector, q = q1 + q2 . The Godambe information
matrix of a information-unbiased composite likelihood is G(θ) = H(θ) = J(θ), and




 Gψψ Gψλ 
,
G(θ) = 


Gλψ Gλλ

where Gψψ is the q1 ×q1 submatrix of G(θ) pertaining to ψ, and Gλλ the q2 ×q2 sub-matrix
of G(θ) pertaining to λ. When λ is unknown, the asymptotic variance of the MCLE of
−1
ψ is given by (Gψψ − Gψλ G−1
λλ Gλψ ) ; when λ is known, the asymptotic variance of the
−1
MCLE of ψ can be shown to be G−1
ψψ . Since Gψλ Gλλ Gλψ is a nonnegative matrix, we have
−1
(Gψψ − Gψλ G−1
≥ G−1
λλ Gλψ )
ψψ . However, the reverse relationship may be observed for an

information-biased composite likelihood, which is illustrated through the equicorrelated
multivariate normal model in the rest of this section. From previous section we know that
the pairwise likelihood CLpair (θ; y) is information-biased for this model.

2.1

Pairwise likelihood in equicorrelated multivariate normal
model

Example 1. Suppose y(1) , . . . , y(n) are n independent observations from the same pdimensional multivariate normal distribution with zero mean and covariance matrix Σ =
σ 2 {(1 − ρ)Ip + ρJp }, where Ip is identity matrix and Jp is a p × p matrix with all entries

6

equal to 1. The common correlation coefficient ρ is the parameter of interest. The
equicorrelated multivariate normal model has been well studied to compare the efficiency
of pairwise likelihood and full likelihood in different settings (Arnold & Strauss, 1991;
Cox & Reid, 2004; Mardia et al., 2009): when σ 2 is unknown, the maximum pairwise
likelihood estimator of ρ, denoted as ρ̂pl , is identical to the MLE of ρ and hence fully
efficient; when σ 2 is known, the maximum pairwise likelihood estimator, denoted as ρ̃pl ,
is less efficient than the maximum likelihood estimator of ρ . Here we are interested in
comparing the asymptotic variances of ρ̂pl and ρ̃pl . The asymptotic variance of ρ̃pl is (
Cox & Reid, 2004)
avar(ρ̃pl ) =

2(1 − ρ)2 c(p, ρ)
,
np(p − 1) (1 + ρ2 )2

(2.1)

where c(p, ρ) = (1−ρ)2 (3ρ2 +p2 ρ2 +1)+pρ(−3ρ3 +8ρ2 −3ρ+2). The asymptotic variance
of ρ̂pl can be shown to be

avar(ρ̂pl ) =

2(1 − ρ)2
{1 + (p − 1)ρ}2 .
np(p − 1)

(2.2)

Comparing the Equations (2.1) and (2.2), we find that as ρ approaches its lower
bound −1/(p − 1), avar(ρ̂pl ) decreases to zero while avar(ρ̃pl ) does not. The ratio of the
asymptotic variances, avar(ρ̃pl )/avar(ρ̂pl ), as a function of ρ is plotted in Figure 1 for
p = 3. We can see that when ρ is positive, ρ̃pl is more efficient than ρ̂pl ; when ρ < 0,
the opposite phenomenon is observed, and when ρ approaches the lower bound −0.5, this
ratio diverges to infinity. We performed the comparisons for different p and observed the
same phenomenon.

7

2.0
1.5
1.0
0.5

The ratio of asymptotic variance

−0.4

−0.2

0.0

0.2

0.4

0.6

0.8

1.0

The value of ρ

Figure 1: The plot of the ratio r(ρ)= avar(ρ̃pl )/avar(ρ̂pl ) at p = 3. The vertical and
horizontal dashed line denotes ρ = 0 and r(ρ) = 1 respectively.
To see that information-biasedness is not a sufficient condition for the paradox to occur,
we consider another information-biased composite likelihood function, the full conditional
likelihood for the same model:

CLF C (θ; y) =

p
Y

f (yr | y(−r) ; θ ),

r=1

where y(−r) denotes the random vector excluding yr . When σ 2 is unknown, the maximum
full conditional likelihood estimator of ρ, ρ̂F C is identical to ρ̂pl and fully efficient (Mardia
et al., 2009); when σ 2 is known, the maximum full conditional likelihood estimator, ρ̃F C is
less efficient than the maximum likelihood estimator for p ≥ 3. Using the formula in Mardia, Hughes, & Taylor (2007), the ratio of the asymptotic variances, avar(ρ̃F C )/avar(ρ̂F C ),
as a function of ρ is plotted in Figure 2 for p = 3. We can see that the ratio is less than
1 for all ρ ∈ [−1/(p − 1), 1].

8

1.0
0.9
0.8
0.7
0.6
0.5
0.4

The ratio of asymptotic variance

−0.5

0.0

0.5

1.0

The value of ρ

Figure 2: The plot of the ratio r(ρ) = avar(ρ̃F C )/avar(ρ̂F C ) at p = 3. The horizontal
dashed line denotes r(ρ) = 1.

2.2

A sufficient condition for the occurrence of paradox

For the inference based on unbiased estimating equations, Henmi & Eguchi (2004) provided a sufficient condition for the occurrence of paradox in the presence of nuisance
parameters, which also applies to the information-biased composite likelihood since its
score function is a special unbiased estimating equation:
Theorem 2.1 Proposition 1.Suppose the composite likelihood function CL(θ | y) is informationbiased and θ = (ψ, λ) with the maximum composite likelihood estimators θ̂c = (ψ̂c , λ̂c ). We
define ψ̃c , the MCLE of ψ when λ is known . Then, ψ̂c has a smaller asymptotic variance
than ψ̃c if λ̂c and ψ̂c are asymptotically independent but λ̂c and ψ̃c are dependent.
2
In the example of equicorrelated multivariate normal model, denote by σ̂pl
the maximum

pairwise likelihood estimator of σ 2 , we can show that the asymptotic covariance between

9

2
is 2ρ(1 − ρ){1 + (p − 1)ρ}σ 2 /(np) which goes to 0 as ρ approaches −1/(p − 1)
ρ̂pl and σ̂pl
2
while the asymptotic covariance between ρ̃pl and σ̂pl
is not equal to zero at ρ = −1/(p−1).

This may explain why the paradox occurs when ρ is close to its lower bound −1/(p − 1).

3

Composite likelihood with more component likelihoods

In this section we consider the inference on a single parameter only, and the simple illustrative examples allow us to calculate the Godambe information matrices or asymptotic
variances analytically.

3.1

Uncorrelated information-biased composite likelihoods

Consider a set of information-unbiased composite likelihood functions CL1 (θ | y), CL2 (θ |
y), . . . , CLm (θ | y) with mutually uncorrelated score functions, it is easy to show that the
product

Qm

i=1

trix G(θ) =

CLi (θ | y) is also information-unbiased and has Godambe information ma-

Pm

i=1

Gi (θ) where Gi (θ) is the information matrix of CLi (θ | y). However, as

shown in the example below, if any composite likelihood CLi (θ | y) is information-biased,
then information additivity may not hold for the product of uncorrelated composite likelihoods.
Example 2. Suppose the random vector (Y1 , Y2 , Y3 ) follows a normal distribution with

10

mean vector µ × (1, 1, 1)T and covariance matrix

 1


Σ=
 ρ


0


ρ
1
0

0 


0 
.


σ2

Assume σ 2 is known, µ and ρ are unknown, and µ is the only parameter of interest.
Suppose we use the independence likelihood function CL12 (µ) = f (y1 ; µ)×f (y2 ; µ), which
is free of the nuisance parameter ρ, to estimate µ. To incorporate the information contained in the independent variable Y3 , we also consider the composite likelihood function
CL123 (µ) = CL12 (µ) × f (y3 ; µ).
It is easy to show that the maximum composite likelihood estimators for CL12 (µ) and
CL123 (µ) are µ̂12 = (ȳ1 + ȳ2 )/2 and µ̂123 = σ 2 (ȳ1 + ȳ2 )/(1 + 2σ 2 ) + ȳ3 /(1 + 2σ 2 ) with
variances (1 + ρ)/2n and [2(1 + ρ)σ 4 + σ 2 ]/n(1 + 2σ 2 )2 respectively. We can compare the
variances of the two maximum composite likelihood estimators directly. For example if
σ 2 = 2, the variance of µ̂123 is (10 + 8ρ)/(25n) which is smaller than (1 + ρ)/(2n) if and
only if ρ > −5/9.
Note that if ρ = −1 this result is expected as (Y1 , Y2 ) determines µ exactly with
µ ≡ (Y1 + Y2 )/2; but the dependence on σ 2 of the range of ρ over which Y3 degrades the
inference is surprising; as σ 2 increases this range approaches [−1, −1/2).

11

3.2

Pairwise likelihood and independence likelihood

Intuitively, a composite likelihood with higher dimensional component likelihoods should
achieve a higher efficiency, although it usually demands more computational cost. In this
subsection we focus on comparing the independence likelihood CLind (θ; y) =
and the pairwise likelihood CLpair (θ; y) =

Qp−1 Qp
r=1

s=r+1

Qp

r=1

f (yr ; θ)

f (yr , ys ; θ). CLpair (θ; y) can be

written as the product of CLind (θ; y) and some pairwise conditional likelihood functions.
Under independence, CLind (θ; y) is identical to the full likelihood, and CLpair (θ; y) =
{CLind (θ; y)}p−1 , which is also fully efficient. For a multivariate normal model with continuous responses, Zhao & Joe (2005) proved that the maximum pairwise likelihood estimator of the regression coefficient has a smaller asymptotic variance than the maximum
independence likelihood estimator.
In fact, within the family of information-unbiased composite likelihood functions pairwise likelihood is at least as efficient as independence likelihood: each bivariate density
f (yr , ys ; θ) has a larger information matrix than f (yr ; θ) and the total number of the bivariate densities in CLpair (θ; y) is p(p−1)/2 > p when p > 2. However, the two composite
likelihoods are information-biased in general especially for complex dependent data and
we may observe the reverse relationship.
A bivariate binary model was used in Arnold & Strauss (1991) to show that the pairwise conditional likelihood could be less efficient than the independence likelihood. For a
bivariate model the pairwise likelihood is the full likelihood and hence fully efficient. Here
we consider a four dimensional binary model which has a complex dependence structure
but also allows us to compare the (asymptotic) variances of different composite likelihood

12

estimators analytically.
Example 3. Suppose (Y1 , Y2 , Y3 , Y4 ) follows a Multinomial(1; θ, θ, θ/k, 1 − 2θ − θ/k),
where k is a positive constant and 0 ≤ θ ≤ k/(2k + 1). The parameter θ controls both the
mean and covariance structures, and we can change the value of k to adjust the strength
of dependence. The value of Y4 is completely determined by 1 − Y1 − Y2 − Y3 . Given a
random sample of size n from this model, we estimate θ based on the independent triplets
(i)

(i)

(i)

(y1 , y2 , y3 ), i = 1, . . . , n.
The full likelihood for the model of (Y1 , Y2 , Y3 ) is

L(θ) =

n
Y

(i)

θ y1

(i)

+y2

i=1

(i)
(i)
(i)
θ
θ (i)
( )y3 (1 − 2θ − )1−y1 −y2 −y3 .
k
k

(3.3)

Solving the score equation we get the maximum likelihood estimator of θ, θ̂ = (ȳ1 + ȳ2 +
ȳ3 )/(2 + 1/k). The exact variance of θ̂ is

Var(

1
θ
ȳ1 + ȳ2 + ȳ3
)= (
− θ2 ).
2 + 1/k
n 2 + 1/k

The independence likelihood function for the model of (Y1 , Y2 , Y3 ) is

CLind (θ) =

n
Y

(i)

(i)

(i)

f (y1 ; θ)f (y2 ; θ)f (y3 ; θ)

i=1

=

n
Y

(i)
(i)
(i)
(i) θ
(i)
(i)
θ
θy1 (1 − θ)y1 θy2 (1 − θ)y2 ( )y3 (1 − )1−y3
k
k
i=1

13

(3.4)

and we can calculate its sensitivity matrix and variability matrix as

Hind (θ) =

2
2
1
1
+
+
+
,
θ 1 − θ kθ k(k − θ)

Jind (θ) =

1
2
2
4
+
−
.
−
2
θ(1 − θ) θ(k − θ) (1 − θ)
(1 − θ)(k − θ)

The pairwise likelihood function is

CLpair (θ) =

n
Y

(i)

(i)

(i)

(i)

(i)

(i)

f (y1 , y2 ; θ)f (y1 , y3 ; θ)f (y2 , y3 ; θ)

i=1

=

n
Y

(i)

θ2y1

(i)

+2y2

(i)

(1 − 2θ)1−y1

(i)

−y2

i=1

(i)
(i)
(i)
θ (i)
θ
( )2y3 (1 − θ − )2−y1 −y2 −2y3 ,
k
k

(3.5)

and we can calculate its sensitivity matrix and variability matrix as

Hpair (θ) =

4
2
4
2(1 + 1/k)2
+
+
+
,
θ kθ 1 − 2θ 1 − (1 + 1/k)θ

θ
θ2
θ
Jpair (θ) = 2A2 θ(1 − θ) + B 2 ( )(1 − ) − 2A2 θ2 − 4AB ,
k
k
k

where A = 2/θ + 2/(1 − 2θ) + (1 + 1/k)/(1 − θ − θ/k), B = 2/θ + 2(1 + 1/k)/(1 − θ − θ/k).
For k = 5, the asymptotic variances of the maximum composite likelihood estimators
for (3.3), (3.4) and (3.5) multiplied by n are plotted as a function of θ in Figure 3.
We can see that when θ < 0.3, the three estimators perform almost equally well; when
θ > 0.3, the full likelihood becomes more efficient than the independence likelihood,
and the independence likelihood estimator is more efficient than the pairwise likelihood
estimator. We also carried out the comparisons for different values of k and found that at
k = 1, both the independence likelihood and the pairwise likelihood are fully efficient, but

14

Figure 3: Left panel: the asymptotic variances (multiplied by n) of the maximum composite likelihood estimators for the full likelihood (solid line), the independence likelihood
(dashed line) and the pairwise likelihood (dotted line). Right panel: the ratio of the
asymptotic variance of pairwise likelihood to the asymptotic variance of independence
likelihood.
when k > 1, the independence likelihood is more efficient than the pairwise likelihood and
the ratio of asymptotic variances approaches 1 when k → ∞. When k < 1, the pairwise
likelihood estimator is more efficient than the independence likelihood estimator and the
ratio of asymptotic variances approaches 1 when k → 0.
This example suggests that in practical applications of composite likelihood inference,
where the models will usually have more complex dependence structure and incomplete
data (e.g., Yi, Zeng, & Cook, 2011), some care is required for the use of higher dimensional
composite likelihood to obtain more efficient estimators. A hybrid composite likelihood
combining lower dimensional marginal and conditional likelihoods with different weights
is suggested to guarantee the improvement of efficiency (Cox & Reid, 2004; Kenne Pagui,
Salvan, & Sartori, 2014a).

15

4

Discussion

As a complement to the discussion on composite likelihood inference in Reid (2012), in
this paper we explored the impact of information bias on the composite likelihood based
inference in different scenarios. An information-unbiased composite likelihood behaves
somewhat like the ordinary likelihood but it can be very inefficient. On the other hand
an information-biased likelihood brings extra difficulty to the computation and is more
likely to exhibit undesirable inferential properties, although the information loss can be
minimized with a set of carefully selected weights.
One way to avoid the paradoxical phenomenon in Section 2 is to convert the composite
score function uc (θ; y) to an unbiased estimating function by projecting (Henmi & Eguchi,
2004; Lindsay, Yi, & Sun, 2011):

u∗c (θ; y) = H(θ)J −1 (θ)uc (θ; y) = arg

min


E ku(θ; y) − ν(θ; y)k2 ,

(4.6)

ν=Auc (θ;y)

where u(θ; y) is the score function of full likelihood, A ranges over all q × q matrices, H(θ)
and J(θ) are the sensitivity matrix and variability matrix. It is easy to check that u∗c (θ; y)
is information-unbiased. Since H(θ) and J(θ) are constant matrices, this projection does
not change the point estimator of θ, and u∗c (θ; y) has the same Godambe information as
uc (θ; y). In the equicorrelated multivariate normal model with θ = (ρ, σ 2 ), the score function of the pairwise likelihood is upl (θ; y) = J(θ)H −1 (θ)u(θ, y) (Pace, Salvan & Sartori,
2011). From Equation (4.6), the projected estimating funtion of upl (θ; y) is equal to the
score function of full likelihood, u(θ; y). In complex models, the required computation for
the projected estimating function u∗c (θ; y) can be intractable and it may be a better idea
16

to design a nuisance-parameter-free composite likelihood function carefully for practical
use. As an example, a pairwise difference likelihood that eliminates nuisance parameters
in a Neyman–Scott problem is described in Hjort & Varin (2008).

References
[1] Arnold, B. C. & Strauss, D. (1991). Pseudolikelihood estimation: some examples.
Sankhyā Series B, 53, 233–243.
[2] Besag, J. E. (1975). Statistical analysis of non-lattice data. The Statistician, 24, 179–
195.
[3] Chandler, R. E. & Bate, S. (2007). Inference for clustered data using the independence
loglikelihood. Biometrika, 94, 167–181.
[4] Cox, D. R. (1975). Partial likelihood. Biometrika, 62, 268–276.
[5] Cox, D. R. & Reid, N. (2004). A note on pseudolikelihood constructed from marginal
densities. Biometrika, 91, 729–737.
[6] Henmi, M. & Eguchi, S. (2004). A paradox concerning nuisance parameters and projected estimating functions. Biometrika, 91, 929–943.
[7] Gao, X. & Song, P. X.-K. (2010). Composite likelihood Bayesian information criteria for model selection in high dimensional data. Journal of the American Statistical
Association, 105, 1531–1540.

17

[8] Godambe, V. P. (1960). An optimum property of regular maximum likelihood estimation. Annals of Statistics, 31, 1208–1212.
[9] Hjort, N. L. & Varin, C. (2007). ML, PL, QL in Markov Chain Models. Scandinavian
Journal of Statistics, 35, 64–82.
[10] Joe, H. & Lee, Y. (2007). On weighting of bivariate margins in pairwise likelihood.
Journal of Multivariate Analysis, 100, 670–685.
[11] Kuk, A. Y. C. (2007). A hybrid pairwise likelihood method. Biometrika, 94, 939–952.
[12] Kenne Pagui, E. C. & Salvan, A. & Sartori, N. (2014a). Combined composite likelihood. To appear in the Canadian Journal of Statistics.
[13] Kenne Pagui, E. C. & Salvan, A. & Sartori, N. (2014b). On full efficiency of the
maximum composite likelihood estimator. Preprint.
[14] Lindsay, B. G. (1982). Conditional score functions:

some optimality results.

Biometrika, 69, 503–512.
[15] Lindsay, B. G. (1988). Composite likelihood methods. Contemporary Mathematics,
80, 221–239.
[16] Lindsay, B. G. & Yi, G. Y. & Sun, J. (2011). Issues and strategies in the selection of
composite likelihoods. Statistica Sinica, 21, 71–105.
[17] Mardia, K. V. & Hughes, G. & Taylor, C. C. (2007). Efficiency of the pseudolikelihood
for multivariate normal and von mises distributions. Preprint.

18

[18] Mardia, K. V. & Kent, J. T. & Hughes, G. & Taylor, C. C. (2009). Maximum likelihood estimation using composite likelihoods for closed exponential families. Biometrika,
96, 975–982.
[19] Pace, L. & Salvan, A. & Sartori, N. (2011). Adjusting composite likelihood ratio
statistics. Statistica Sinica, 21, 129–148.
[20] Reid, N. (2012). Likelihood inference for complex settings. Canadian Journal of
Statistics, 40, 731–744.
[21] Varin, C. (2008). On composite marginal likelihoods. Advances in Statistical Analysis,
92, 1–28.
[22] Varin, C. & Reid, N. & Firth, D. (2011). An overview of composite likelihood methods. Statistica Sinica, 21, 5–42.
[23] Varin, C. & Vidoni, P. (2005). A note on composite likelihood inference and model
selection. Biometrika, 92, 519–528.
[24] Xu, X. & Reid, N. (2011). On the robustness of maximum composite likelihood
estimate. Journal of Statistical Planning and Inference , 141, 3047–3054.
[25] Yi, G. Y. & Zeng, L. L. & Cook, J. R. (2011). A robust pairwise likelihood method for
incomplete longitudinal binary data arising in clusters. Canadian Journal of Statistics,
39, 34–51.
[26] Zhao, Y. & Joe, H. (2005). Composite likelihood estimation in multivariate data
analysis. Canadian Journal of Statistics, 33, 335–356.

19

