J. R. Statist. Soc. B (2017)

A minimum relative entropy based correlation model
between the response and covariates
Bhaskar Bhattacharya
Southern Illinois University, Carbondale, USA

and Mohammad Al-talib
Yarmouk University, Irbid, Jordan[Received May 2015. Final revision June 2016]
Summary. A semiparametric model is presented utilizing dependence between a response
and several covariates. We show that this model is optimum when the marginal distributions of
the response and the covariates are known. This model extends the generalized linear modeland the proportional likelihood ratio model when the marginal distributions are unknown. New
interpretations of known models such as the logistic regression model, density ratio model andselection bias model are obtained in terms of dependence between variables. For estimation ofparameters, a simple algorithm is presented which is guaranteed to converge. It is also the same
regardless of the choice of the distribution for response and covariates; hence, it can fit a very
wide variety of useful models. Asymptotic properties of the estimators of model parameters arederived. Real data examples are discussed to illustrate our approach and simulation experiments
are performed to compare with existing procedures.
Keywords: Dependence; Information; Logistic regression; Semiparametric generalized linear
model

1.

Introduction

Luo and Tsai (2012) described neuropsychological scale data where the response variable is thescore from the trail making test (part A) measuring 334 patients’ processing speed in seconds,
and the covariates are years of education, age and diagnosis. Often these types of data have several covariates (Figs 1 and 2 in Section 7 show scatter plots of scores versus years of education
and age), have unknown statistical distributions and known statistical procedures fail to work
properly. Stamey et al. (1989) examined the correlations between the level of prostate-speciﬁc
antigen and several clinical measures in 97 men who were about to receive a radical prostatectomy. The goal (Hastie et al., 2009) is to predict the logarithm of prostate-speciﬁc antigen level,lpsa, from a number of measurements including log-cancer-volume, lcavol, log-prostate-weight,
lcp, age and logarithm of capsular penetration, lcp (Figs 5–7 in Section 7 show scatter plotsof lpsa versus lcavol, lweight and lcp respectively). Although there are moderately high correlations between the response and most covariates, a linear regression model does not consider
the non-linearity at the edges of the data. In both examples, our goal is to develop a modelfor the response variable by utilizing its dependence on the covariates while maintaining all the

Address for correspondence: Bhaskar Bhattacharya, Department of Mathematics, Southern Illinois University,
Carbondale, IL 62901-4408, USA.
E-mail: bhaskar@siu.edu© 2016 Royal Statistical Society

1369–7412/17/79000

2

B. Bhattacharya and M. Al-talib

marginal information about the response and the covariates that is obtained from the sample
in the model.
For a response Y with covariates Zi , 1  i  d, the correlation coefﬁcient corr.Y , Zi / = cimeasures any linear relationship between Y and Zi . When |ci | is much smaller than 1, thenthe linear relationship between Y and Zi is very weak. Then we may think of the relationshipbetween Y and Zi s as non-linear in nature. In this sense, ci not only measures the strengthof the linear but also the non-linear relationship between Y and Zi with higher or lower ci
referring respectively to more linearity or more non-linearity. Thus ci can be considered as a
global measure of dependence between Y and Zi . For example, if E.Y |Zi / = Zi2 with Zi being
symmetrically distributed around 0, then ci = 0, then the dependence of Y on Zi is entirely
non-linear. Fig. 8(c) in Section 8 considers a similar situation with E.Y |Z/ = |Z|.
In this paper, we develop a model for Y based on corr.Y , Zi / = ci , 1  i  d. We shall see
that, when the speciﬁed ci s are high, the model that is obtained is almost linear (e.g. prostatedata), whereas, when the speciﬁed ci s are low, the model that is obtained is non-linear (e.g.
the trail making data). In any scientiﬁc study, the initial choice of covariates from a vast pool
may be difﬁcult. Given a set of covariates, transformations of covariates or increasing the number of correlation constraints with original or transformed variables might prove useful for abetter ﬁt of the model developed. See Sections 7 and 8 for discussions on transformations ofcovariates.
When the marginal distributions of Y and the Zi s are known, our procedure has close connections with the maximum entropy (ME) principle, which may be stated as follows:‘when selecting a model for a given situation it is often appropriate to express the prior information interms of constraints. However, one must be careful so that no information other than these speciﬁed
constraints is used in model selection. That is, other than the constraints that we have, the uncertainty
associated with the probability distribution to be selected should be kept at its maximum’ (Jaynes,
1957).

The ME principle can be generalized to the concepts of Kullback–Leibler (KL) distance andI-projection, as deﬁned below (Csiszár, 1975). For two probability measures Q and P, the KL
distance (or, relative entropy) between Q and P is deﬁned as   dQ ln
dQ,
if Q  P,
I.Q|P/ =
.1:1/
dP
∞,
otherwise:(Q  P means that Q is absolutely continuous with respect to P.) Although I.Q|P/ is not a
metric, it is always non-negative and equals 0 if and only if Q = P. Hence it is often interpreted
as a measure of ‘divergence’ or ‘distance’ between Q and P. For a given P and a speciﬁed set ofprobability measures C, it is often of interest to ﬁnd the QÅ ∈ C which satisﬁes
I.QÅ |P/ = inf I.Q|P/Q∈C.< ∞/:.1:2/

Such a QÅ is called the I-projection of P onto C. Csiszár (1975) has shown that QÅ exists
uniquely if C is convex and variation closed and there is a Q ∈ C such that I.Q|P/ < ∞ (when P
is uniform, equation (1.2) becomes the ME principle).
Associations between random variables have been of interest to statisticians and probabilists
over several recent decades (Agresti, 2013). Recently, Reshef et al. (2011) have proposed ameasure of dependence for two-variable relationships, known as the maximal information coefﬁcient. Although the maximal information coefﬁcient captures a wide range of associations
both functional and not, it is restricted to only two variables at a time. Szekely et al. (2007,

Minimum Relative Entropy Correlation Model32009) have proposed distance correlation measures of dependence between variables. But thispaper proposes models utilizing dependence based on correlations between variables and their
marginal distributions.
In Section 2 ﬁrst we show that the model proposed is optimum when the marginal distributions of Y and Z are known. Then we propose a model for Y when the marginal distributionsof Y and the Zi s are unknown, and we discuss its properties. In Section 3, we describe the relationship of the model with other existing models. Also, new interpretations of known models
are obtained. In Section 4, we present an algorithm to estimate the parameters of the model
proposed. This algorithm is guaranteed to converge and is applicable to any distributions of Y
and Zi s. Consistency and asymptotic normality of the parameters are addressed in Section 5.
In Section 6, we present results from two simulation experiments, for discrete and continuous
response cases. In Section 7, two real data examples are considered: one from neuropsychological test scores, and the other involving prostate cancer trials. Discussion and ﬁnal comments
are in Section 8. Further details on duality, proofs of theorems and residual plots are in the
supplemental ﬁle that is available on line.
The data that are analysed in the paper and the programs that were used to analyse them canbe obtained fromhttp://wileyonlinelibrary.com/journal/rss-datasets
2.

Model proposed and its properties

Let the response variable Y have cumulative distribution function (CDF) P1 .y/, and marginal
probability density functions (PDFs) p1 .y/. Those for a d-dimensional covariate vector Z =.Z1 , : : : , Zd / are P2 .z/ and p2 .z/ respectively. First we assume that the marginal distributions of
Y and Z are known. As we are interested in the dependence between Y and Z, let .Y , Z/ have a(unknown) joint CDF Q.y, z/. Let the marginal CDFs of Y and Z that are derived from Q.y, z/
be Q1 .y/ and Q2 .z/ (PDFs q1 .y/ and q2 .z/) respectively. Now consider a class C of all joint
probability distributions Q.y, z/ of .Y , Z/ whereC = {Q : corr.Y , Zi / = ci , 1  i  d, Q1 .y/ = P1 .y/, ∀ y ∈ R, Q2 .z/ = P2 .z/, ∀ z ∈ Rd },.2:1/

for a given vector of constants c = .c1 , : : : , cd /, − 1  ci  1. Using the marginal distributions of
Y and Z, we can equivalently express
C = {Q : E.Y Z/ = c , Q1 .y/ = P1 .y/, ∀ y ∈ R, Q2 .z/ = P2 .z/, ∀ z ∈ Rd },.2:2/√
where c = .c1 , : : : , cd /, ci = ci {var.Y/ var.Zi /} + E.Y/ E.Zi /, ∀ i. Without loss of generality,
we shall use ci = ci .
The independence model P = P1 P2 ignores any dependence relationship between Y and Z.
We propose a model QÅ for .Y , Z/ to be that Q in C that is closest to P in the minimum KL
divergence (recall the ME principle). Theorem 1 below shows that the criteria in C lead us to
the solution given by (where T denotes transpose)qÅ .y, z/ = 

exp{βÅT zy + hÅ1 .y/ + hÅ2 .z/} p1 .y/ p2 .z/
exp{βÅT vu + hÅ1 .u/ + hÅ2 .v/} p1 .u/ p2 .v/dudv

:

.2:3/

Here the (adjustment) functions hÅ1 .y/ and hÅ2 .z/ are needed so that the solution maintains theY - and Z-marginals as speciﬁed in C. Sometimes the exact expressions of hÅ1 .y/ and hÅ2 .z/ areavailable, as in example 1 below. In other cases they are approximated by using the algorithm

4B. Bhattacharya and M. Al-talib

given. Deriving the conditional distribution from the joint distribution (2.3), a model of Y , for
given Z, isqÅ .y|z/ = 

exp{βÅT zy + hÅ1 .y/} p1 .y/
exp{βÅT zu + h1 .u/} p1 .u/du

:

.2:4/

Since theorem 1 characterizes equation (2.3) as an optimum solution under C, we shall saythat equation (2.4) is the optimum model for Y using Z (meaning, derived from equation (2.3)).
The following theorem proves the optimality of equation (2.4) by using the duality techniquefrom Bhattacharya and Dykstra (1995, 1997). More details on duality are in the supplemental
ﬁle that is available on line.
Theorem 1. Assume that c, p1 .y/ and p2 .z/ in C are known. Then model (2.4) is the optimum
model of Y for given Z that incorporates the correlation and marginal information in C, in thesense that the corresponding joint distribution QÅ (in C) is the closest to the independence
model P in the minimum KL distance than any other Q in C.
The following two examples use theorem 1 to prove a characterizing property of the multivariate normal distribution in the sense that the closest distribution from the product of (multivariate) normal distribution marginals onto the set C is also multivariate normal. Example 1
shows exact expressions for hÅ1 .y/ and hÅ2 .z/. Comparing these two examples, we can see theimportance of the presence of marginal constraints.
2.1. Example 1
For simplicity, we consider the case d = 2, and one can extend it to any d. Consider the
random vector .Y , Z1 , Z2 / such that the marginals are Y ∼ P1 = N.0, 1/, and .Z1 , Z2 / ∼ P2 =N 2 {0, .1, ρ==ρ, 1/}, − 1  ρ  1, for a known value of ρ. Consider the class C1 of probability
distributions Q1 where
C1 = {Q1 : EQ1 .Y Z1 / = ρ1 , EQ1 .Y Z2 / = ρ2 , Y ∼ N.0, 1/, .Z1 , Z2 / ∼ N {0, .1, ρ==ρ, 1/}},and −1  ρ1 , ρ2  1 are known. Assume that the values of ρ, ρ1 and ρ2 are such that Σ1 derived
below is positive deﬁnite (as the numerical example shows). We look for QÅ1 in C1 which is
closest to P = P1 P2 in the minimum KL divergence (P1 and P2 independent). From equation(2.3) and theorem 1, the solution is QÅ1 , where q1Å .y, z/ = dQÅ1 =dμ (μ is Lebesgue measure), is
given by (z = .z1 , z2 / and v = .v1 , v2 //Å yz + β Å yz + hÅ .y/ + hÅ .z/} p .y/ p .z/
exp{β11
1
1
2
12 2
1
2
:
.2:5/
q1Å .y, z/ = ÅÅÅÅ
exp{β11 uv1 + β12 uv2 + h1 .u/ + h2 .v/} p1 .u/ p2 .v/dudv
Consider the case when .Y , Z1 , Z2 / ∼ N3 .0, Σ1 /, where Σ1 = .σij /, with σii = 1, σ12 = ρ1 , σ13 =ij
ρ2 , σ23 = ρ and σij = σji , ∀ i, j. Then from the expression for Σ−1
1 = .σ /, the joint PDF of.Y , Z1 , Z2 / isq.y, z1 , z2 / =exp{−1=.2k1 /}1=2.2π/3=2 k1{.1 − ρ2 /y2 + .1 − ρ22 /z21 + .1 − ρ21 /z22 − 2.ρρ2 − ρ1 /yz1

−2.ρρ1 − ρ2 /yz2 − 2.ρ1 ρ2 − ρ/z1 z2 },where k1 = det.Σ1 / = 1 − ρ2 − ρ21 − ρ22 + 2ρρ1 ρ2 .

.2:6/

Minimum Relative Entropy Correlation Model

5

The distribution N3 .0, Σ1 / is in C1 . Using the normal PDF expressions for p1 .y/ and p2 .z/
in equation (2.5) and comparing the resulting exponent with that in equation (2.6), we can see
using algebra that these two exponents would be the same (except for the additive constants) if
we set−ρ21 − ρ22 + 2ρρ1 ρ2 2.y − 1/,2k1−ρ2 − ρ21 + 2ρρ1 ρ2 2
−ρ2 − ρ22 + 2ρρ1 ρ2 2
.z1 − 1/ +.z2 − 1/h2Å .z1 , z2 / =2
2k1 .1 − ρ /2k1 .1 − ρ2 /−2ρ + 2ρ3 + ρ1 ρ2 + ρρ21 + ρρ22 − 3ρ2 ρ1 ρ2+.z1 z2 − ρ/
k1 .1 − ρ2 /
Å = .ρρ − ρ /=k and β Å = .ρρ − ρ /=k . Also,  hÅ .y/ dP .y/ =  hÅ .z/ dP .z/ = 0 isand β11
2
1
1
1
2
1
1
2
12
1
2
satisﬁed (see the proof of theorem 1). Since the solution of equation (1.2) is unique, QÅ1 =N3 .0, Σ1 /[∈ C] must be the solution to our problem. As a numerical example, let ρ = 0:7, ρ1 =
Å = −0:132 and β Å = 0:763.
0:3 and ρ2 = 0:5; then q1Å has β1112h1Å .y/ =2.2. Example 2
Consider the same setting as in example 1. Removing the marginal information from C1 , consider
the class C2 of probability distributions Q2 where
C2 = {Q2 : EQ2 .Y Z1 / = ρ1 , EQ2 .Y Z2 / = ρ2 },where −1  ρ1 , ρ2  1 are known. We look for QÅ2 in C2 which is closest to P = P1 P2 (with thesame P as in example 1) in the minimum KL divergence. Using a similar technique to those intheorem 1 and its proof, the solution is Q2Å , where q2Å .y, z/ = dQ2Å =dμ (μ is Lebesgue measure)is given by
q2Å .y, z/ = =Å yz + β Å yz / p .y/ p .z , z /exp.β21
1
1
2 1 2
22 2
Å uv + β Å uv / p .u/ p .v , v /dudv
exp.β21
1
1
2 1 2
22 2
Å yz + β Å yz − y2 =2 − [1={2.1 − ρ2 /}].z2 − 2ρz z + z2 //exp.β21
1
1 2
1
2
22 2
Å uv + β Å uv − u2 =2 − [1={2.1 − ρ2 /}].v2 − 2ρv v + v2 //dudv
exp.β21
1
1 2
1
2
22 2

:

.2:7/

Write the exponent of the numerator of the rightmost term of equation (2.7) as − 21 .y, z1 ,T
z2 /Σ−1
2 .y, z1 , z2 / , where
⎞
⎛
Å
Å
−β22
1
−β21
⎝ −β Å
Σ−1
1=.1 − ρ2 / −ρ=.1 − ρ2 / ⎠,2 =21
Å−β22 −ρ=.1 − ρ2 /1.1 − ρ2 /⎛⎞ÅÅ ρ + βÅÅβ21
1
β21 + β22 ρ
22
1 ⎝ Å
Åρ
Å
Å βÅ ⎠
Σ2 =β21 + β22
1 − .1 − ρ2 /β22
ρ + .1 − ρ2 /β21
22
k2
Å
Å
Å
Å
Å
2
2
β21 ρ + β22 ρ + .1 − ρ /β21 β22
1 − .1 − ρ /β21
Å − β Å − 2β Å β Å ρ. Setting the .1, 2/th and .1, 3/th elements of Σ equal to
and k2 = 1 − β21
2
22
21 22
Å = k .ρ − ρρ /=.1 − ρ2 / and β Å = k .ρ − ρρ /=.1 − ρ2 /.
ρ1 and ρ2 respectively, we obtain β21
2 1
2
2 2
1
22

6

B. Bhattacharya and M. Al-talib

Å and β Å in the expression for k above, we obtain the quadratic k2 .ρ2 +Using these values of β21
2
2 1
22
2
ρ2 − 2ρρ1 ρ2 / + k2 − 1 = 0, which is easily solved for k2 .
Since Q2Å = N.0, Σ2 / ∈ C2 , and the solution of equation (1.2) is unique, QÅ2 must be thesolution to our problem. Note from Σ2 above that the marginal distributions of Y and Z are
now different from those of the solution in example 1 (Σ1 ). As a numerical example, using
Å = −0:076, β Å = 0:442 andthe same ρ, ρ1 and ρ2 as in example 1, we obtain that q2Å has β2122k2 = 0:777 (discarding the negative value).
Å , β Å / = .β Å , β Å /. Thus, ﬁxing the marginals
Comparing these two examples, C1 ⊂ C2 and .β1112
21 22
in equation (2.2) makes a difference in the solution of equation (1.2). It may seem difﬁcult tomotivate equation (2.7) with regard to P = P1 P2 when the marginal distributions are not speciﬁedin C2 ; however, as shown in Section 3.2, this is the scenario that produces the proportional
likelihood ratio model.
Analytical solutions for hÅ1 .y/ and hÅ2 .z/ can be obtained by using a general algorithm given in
Bhattacharya (2006). For a given sample, the variables are discrete, and we provide an algorithm
to ﬁnd h1Å .y/ and hÅ2 .z/ in this paper. In constructing a model, we would beneﬁt from utilizing
as much dependence between Y and Z as possible. So, if corr(Y , Z) is low, then we could
employ some functions, e.g. g1 .Y/ and g2 .Z/, provided that corr{g1 .Y/, g2 .Z/} is higher. Similar
approaches have been used by Fokianos and Kaimi (2006) and Guerrero and Johnson (1982).
We discuss this further in Sections 3, 7 and 8.
For an unknown data-generating process, however, the true values of c, p1 .y/ and p2 .z/ will
not be known. Then we can neither deﬁne C as in expression (2.1); nor does the optimalityof equation (2.4) hold as above. However, for large sample size, the marginal distributions of
Y and Z can be estimated by using their empirical distributions along with the population
correlations between Y and Zi by the corresponding sample values. Using these observed valuesof constraints from the sample, we deﬁne a discrete version of C (namely, K in expression (4.2)) in
the same way as in equation (2.1) in Section 4. Assuming that these sample (moment) constraints
represent the corresponding population (moment) constraints efﬁciently and consistently, theresulting model (though not necessarily optimum) is expected to perform well. Model (2.4) now
becomes semiparametric because, along with the parametric component β, there is also thenon-parametric component p1 .y/, where no distributional assumption about p1 .y/ is made.
We would like to begin with as much information as possible about the response and covariatesthat is available but with independence between them, i.e. ﬁrst ﬁnd P̂ = P̂ 1 P̂ 2 . This choice of P̂
also gives us all known models as shown in Section 3. If P̂ were to be a uniform distribution(ME formulation), that would mean throwing away some of the observed information about
the response and covariates.
Next we discuss some properties of equation (2.4) with p1 .y/ unknown. The likelihood ratioisqÅ .y2 |z/
p1 .y2 /log Å= βT z.y2 − y1 / + h1 .y2 / − h1 .y1 / + log,.2:8/
q .y1 |z/p1 .y1 /for any z, y1 and y2 such that p1 .y1 / > 0 and p1 .y2 / > 0. If z is one dimensional, thenβ = logq.y + 1|z + 1/=q.y|z + 1/,
q.y + 1|z/=q.y|z/i.e. exp.β/ measures the likelihood ratio that the response increases by 1 unit when the covariate
increases by 1 unit. Hence equation (2.4) can be called a (generalized) proportional likelihood
ratio model as the likelihood ratio (2.8) depends on z through βT z (we call it ‘generalized’
because of the extra term h1 .y2 / − h1 .y1 /, compared with Luo and Tsai (2012)).

Minimum Relative Entropy Correlation Model7

Let θ = βT z. If μ.θ/ = E.Y |Z = z/, then var.Y |Z = z/ = μ .θ/, which is a property of thesufﬁcient statistics in exponential families (which model (2.4) is not when p1 .y/ is unknown).
Also, it can be easily seen that β is invariant to monotone increasing transformations on Y andlocation shifts of Z.
3.

Relationship to other models

3.1. Logistic regression model
We show below that model (2.3) reduces to the logistic regression model when Y is binary. To
see this, let Y take values y1 = 0 or y2 = 1 and Z be ﬁxed at a value zj , where j = .j1 , : : : , jd /. Also,
let {rij } be given constants, where Σi,j rij = 1. Model (2.4) is given by
rij exp.βT yi zj + βk+1,i /P.Y = yi |Z = zj / =2
w=1

rwj exp.βT yw zj + βk+1,w /so deﬁning π.zj / = P.Y = 1|Z = zj / we obtain
π.zj / =r2j exp.βT zj + βk+1,2 /r1j exp.βk+1,1 / + r2j exp.βT zj + βk+1,2 /

:

Then,
π.zj /
r2j exp.βk+1,2 − βk+1,1 /exp.βT zj /=1 − π.zj /r1jsologit{π.zj /} = β0 + βT zj ,which is the logistic regression model, where 
r2j
P.Y = 1/
+ βk+1,2 − βk+1,1 = ln
β0 = lnr1j
P.Y = 0/+ βk+1,2 − βk+1,1

depends on the (unknown) marginal distribution of Y . By theorem 1, if ci s and the marginal
distributions of Y and Z as in expression (2.2) were speciﬁed, then the logistic regression model(is the conditional distribution of Y for given Z = z of a joint distribution in C, which) is theclosest to the independence model subject to the criteria in C.
3.2. Proportional likelihood ratio model
The semiparametric proportional likelihood ratio model that was proposed by Luo and Tsai(2012) (see also Chan (2013)) is given byÅ .y|z/ = qPLR

exp.βT zy/p1 .y/
exp.βT zu/p1 .u/du

,

.3:1/

where the Y -marginal p1 .y/ is not speciﬁed. Note that equation (3.1) does not carry the term
h1 .y/ in the exponent like in equation (2.4). Hence the marginal PDF of Y under model (3.1) is
exp.βT zy/ p1 .y/p2 .z/dz
,T
exp.β zu/p1 .u/du8

B. Bhattacharya and M. Al-talib

which is not necessarily equal to p1 .y/. The estimation of parameters of model (3.1) depends
on an algorithm; however, the right condition to guarantee the convergence of the algorithm is
not tractable. Also, the covariates Zi need to be bounded.
In fact, following Section 2 and as in example 2, the I -projection of P onto C = {Q :corr.Y , Zi / = ci , 1  i  d} is QÅ , wheredQÅ= qÅ .y, z/ = dμexp.βÅT zy/p1 .y/p2 .z/
exp.βÅT vu/ p1 .u/ p2 .v/dudv,

from which equation (3.1) follows. Thus if p1 .y/ and p2 .z/ were known, then model (3.1) is optimum in the KL sense subject to C . Leaving p1 .y/ and p2 .z/ unknown, model (3.1) is a special
case of model (2.4) setting hÅ1 .y/ = hÅ2 .Z/ = 0. Note that there can be many candidate distributions of Y and Z that can produce the same level of dependence c. It is shown in examples 1 and2 that absence of the marginals in the constraints may lead to a different model. It is our pointthat by setting the restrictions as C instead of C we would capture the marginal informationof Y and Z in the model, which might provide more accurate inference about Y by using Z.
3.3. Generalized linear model
Consider a generalized linear model that speciﬁes the conditional distribution of Y given Z as
fglm .y|z/ = exp

θy − b.θ/+ c.y, τ / ,a.τ /where a, b and c are given functions so that fglm .y|z/ is a PDF. Now let θ = βT z and set c† .y, τ / =h1 .y/ + log{p1 .y/}, β† = β=a.τ / and†Tb .θ/ = a.τ / ln exp{β zu + h1 .u/}p1 .u/ du :Then†fglm .y|z/ = exp β†T zy −b† .θ/+ c† .y, τ /
a.τ /

†

is of the form (2.4). Thus, fglm .y|z/ extends the generalized linear model by allowing the partp1 .y/ in c† .y, τ / to be unspeciﬁed. Since p1 .y/ is unspeciﬁed, the standard methods of statisticalinference for exponential families do not apply.
3.4. Density ratio models
In the presence of case–control data, the two-sample density ratio model (Qin, 1998) is f.z|y =1/ = exp{α + βT h.z/}f.z|y = 0/, for some functions h. Using the Bayes formula, we can write
this asexp{βT h.z/y}p1 .y/
qÅ .y|z/ = :
exp{βT h.z/u}p1 .u/duThus using theorem 1 when the distributions of Y and Z are known, we can characterize it as
the optimum in the sense that it is the conditional distribution derived from a joint distribution
in the set C = {Q : corr{Y , h.Z/} = c} which is closest (KL) to the independence model P. Whenthe distributions of Y and Z are unknown, it is a special case of model (2.4) with no marginal

Minimum Relative Entropy Correlation Model

9

information preserved. Kay and Little (1987) and Guerrero and Johnson (1982) have consideredtransformations of the explanatory variables in logistic regression models for binary data.
Also, the m-sample density ratio model of Fokianos and Kaimi (2006) can be characterized
as the best in the sense that it is the conditional distribution derived from a joint distribution in
the set C = {Q : corr{ki .Y/, hj .Z/} = cij , ∀ i, j} which is closest (KL) to the independence model,
assuming that the distributions of Y and Z are known. When they are unknown, the m-sample
density ratio model is a special case of model (2.4) with no marginal information preserved.
3.5. Dependence measuresJoe (1989) and Good (1976) suggested the use of δX1 ,:::,Xm = I.fX1 ,:::, Xm |Πmi=1fXj / as a measure ofmultivariate dependence for random vectors X1 , : : : , Xm . Joe (1989) showed that, for multivariate
normal distributions with m = 2, δX1 ,X2 is equivalent to |ρ|, where ρ is the correlation coefﬁcient.
Thus, in general, δX1 ,:::,Xm can be useful as a measure of functional dependence between random
vectors X1 , : : : , Xm . In this paper we identify the joint distribution fY ,Z1 ,:::,Zd that is closest to
the independence model fY fZ1 ,:::,Zd by minimizing I.fY ,Z1 ,:::,Zd |fY fZ1 ,:::,Zd / under dependence
and marginality restrictions.
Interpretations similar to those in Sections 3.1–3.4 can also be obtained for other models suchas the Gilbert–Lele–Vardi biased sampling model (Gilbert et al., 1999) and the semiparametric
single-index model (Ichimura, 1993). In this context, note that Rathouz and Gao (2009) and
Huang and Rathouz (2012) modelled the mean directly, where the parameter β is essentially acontrast in the mean response.
4.

Estimation

Suppose that a random sample .Yi , Zi /, i = 1, : : : , n, is available from Q, the joint distributionof .Y , Z/. Let y1 , : : : , yk be the distinct observed values of Y , where yi occurs ni times, 1  i  k,
Σki=1 ni = n. Similarly, for Z = .Z1 , : : : , Zd /, let {.zs1 , : : : , zsks /} be the distinct observed values of
Zs , s = 1, : : : , d. The index j = .j1 , : : : , jd / runs lexicographically from .1, : : : , 1/ to .k1 , : : : , kd /.
Also, let .z1j1 , : : : , zdjd / occur mj times and, jointly, .yi , z1j1 , : : : , zdjd / occur nij times, 1  i  k,
1  jt  kt , 1  t  d, Σj mj = Σi,j nij = n.
For discrete Y and Z, we write Q as q = .qij / and P as p = .pij /. The observed CDFs areP̂ = .p̂ij / and Q̂ = .q̂ij /. Here qij = P.Y = yi , Z = zj / = q.yi , zj /. Then the observed joint probability
mass function is q̂ij = P.Y = yi , Z1 = z1j1 , : : : , Zd = zdjd / = nij =n, which is used to calculate theobserved value of cs = Σij yi zsjs nij =n, s = 1, : : : , d, and the marginal distributions of Y and Z asΣj qij = ni =n, i = 1, : : : , k, Σi qij = mj =n, j = .j1 : : : jd / = .1, : : : , 1/, : : : , .k1 : : : kd / respectively.
When the CDFs P and Q in equation (1.2) are replaced by the corresponding empirical CDFs
P̂ = .p̂ij / and Q̂ = .q̂ij / respectively, the problem reduces to the discrete case and can be expressed
as
 
q̂ij
inf I.Q̂|P̂/ = inf,.4:1/
q̂ij lnp̂ij.q̂ij /∈K ij
Q̂∈K
where the constraints in C in expression (2.2) are replaced by their discrete versions in K inequation (4.1), which is deﬁned asmjni
K = q = .qij / : yi zsjs qij = cs , s = 1, : : : , d, qij = , i = 1, : : : , k, qij = ,
n
n
i
ij
j

.4:2/
j = .j1 : : : jd / = .1, : : : , 1/, : : : , .k1 : : : kd / :10

B. Bhattacharya and M. Al-talib

For easier exposition, we number k1 : : : kd of j-values using t = 1, : : : , k1 : : : kd (as in expression(4.3) below), so t = j for some j, by slight abuse of notation. In equation (4.1), P̂ = .p̂ij / is the
sample estimate of P = P1 P2 , so p̂ij = .ni =n/mj =n, and the estimate to be found is Q̂ = .q̂ij /.
Although P̂ satisﬁes the row and column constraints, it may not satisfy constraints Σij yi zsjs qij= cs . For solving problems such as equation (4.1) for multiway contingency tables, Bhattacharya and Dykstra (1997) have shown that a Fenchel duality theorem can be used to identify
a dual optimization problem, which might be easier to solve. When the constraint region is anintersection of convex cones, they proposed an algorithm which is guaranteed to converge. To
adopt this strategy in our case, for a given n, consider problem (4.1) as the primal problem.
For simplicity of notation, we suppress the dependence on n (for sets and dual solutions hij )up to expression (4.10) below.
First we express equation (4.2) as an intersection of convex sets as follows:
k :::k1 d
K = .∩ds=1 Ks / ∩ .∩kr=1 Mr / ∩ .∩t=1
Nt /,

where



Ks = q = .qij / : .yi zsjs − cs /qij = 0 ,ijMr = q = .qij / : .Iijr − nr =n/qij = 0 ,ijNt = q = .qij / : .Jijt − mt =n/qij = 0 ,1  s  d,r = 1, : : : , k,

⎫
⎪
⎪
⎪
⎪
⎪
⎬

⎪
⎪
⎪
⎪
⎭
t = 1, : : : , k1 : : : kd , ⎪

.4:3/

ij

where Iijr = 1 if r = i, and Iijr = 0 otherwise, Jijt = 1 if t = j, and Jijt = 0 otherwise. These will bereferred to as K-constraints, M-constraints and N-constraints respectively. The dual cone KÅof K is given by (Rockafellar, 1970)
k :::kKÅ = .⊕d KÅ / ⊕ .⊕k MÅ / ⊕ .⊕ 1 d NÅ /,
.4:4/
s

s=1

r=1

r

t=1

t

where ‘⊕’ indicates that the direct sum of vectors KÅs , MÅr and NÅt are the dual cones of
Ks , Mr and Nt respectively, each consisting of vectors of length kk1 : : : kd :KsÅ = h1s = .h1sij / = .γ1s .yi zsjs − cs //, γ1s ∈ R, ∀ i, j ,MrÅ = h2r = .h2rij / = .γ2r .Iijr − nr =n//, γ2r ∈ R, ∀ i, j ,NÅ = h = .h / = .γ .J t − m =n//, γ ∈ R, ∀ i, j :
t

3t

3t

3tij

ij

3t

t

Then h = .h / ∈ KÅ is a vector of length kk : : : k , where each of its elements is a sum of three
1

ij

d

components, and can be written as (here we suppress n and v as hij = hnvij , h1sij = hnv1sij , etc.)
hij =

d

h1sij +

s=1

k
r=1

h2rij +

k1:::kd

and the dual problem to problem (4.1) is given byinfÅ p̂ij exp.hij /:h∈K ijh3tij ,

.4:5/

t=1

.4:6/

For ﬁxed n, problem (4.1) is solved by iteratively ﬁnding the I -projections onto the K-, Mand N-constraints. Each cycle has d + k + k1 : : : kd steps, which are the same as the numberof constraints. First for the K-constraints, for 1  s  d, if, at the vth cycle, the .s − 1/th step(referred to as the .v, s − 1/th step), the solution is q̂.n, v, s − 1/ = .q̂ij .n, v, s − 1//, then, at the.v, s/th step, we solve for (q = q̂.n, v, s/) in.4:7/
inf I{q|q̂.n, v, s − 1/}:q∈Ks

Minimum Relative Entropy Correlation Model

11

The dual objective function at the .v, s/th step is obtained from the convex conjugate (see theon-line supplement, theorem S1) of the function in expression (4.7). The corresponding dual
problem solves
inf

h1sij ∈KÅs ij

q̂ij .n, v, s − 1/ exp.hij / = inf

γ1s ∈R ij

q̂ij .n, v, s − 1/ exp.hij /:.4:8/

Note that the dual problem (4.8) amounts to solving for the scalar γ1s only, which is considerablyeasier than solving for q in expression (4.7). The solution to the right-hand side of equation (4.8),
say, at γ1s = γ1nvs is obtained by differentiating it with respect to γ1s and setting it equal to 0,i.e.
q̂ij .n, v, s − 1/.yi zsjs − cs / exp.hij / = 0,
.4:9/
ij

which can be solved easily by the Newton-Raphson method. Then, using theorem S1, we obtainthe solution to expression (4.7) asq̂ij .n, v, s − 1/ exp.h1sij /q̂ij .n, v, s/ =,.4:10/
q̂ab .n, v, s − 1/ exp.hab /abwhere b = .b1 , : : : , bd /.
Next, for the M-constraints, we must match the y-marginals of the current table to the observed values nr =n. If, at the .v, d + r − 1/th step, the solution is q̂.n, v, d + r − 1/, then theexact dual solution at the .v, d + r/th step is γ2nvr = log[nr ={nq̂r−1,+ .n, v, d + r − 1/}], whereq̂r−1,+ .n, v, d + r − 1/ = Σj q̂r−1,j .n, v, d + r − 1/. This amounts to rescaling the Y -marginals of
the current table as in step 3 of the algorithm (see below).
Finally, considering the N-constraints, we must match the z-marginals of the current table to the observed values mj =n; the exact dual solution at the .v, d + k + t/th step is γ3nvt =log[mj ={nq̂+,t−1 .n, v, d + k + t − 1/}], (recall that t = j for some j) where q̂+,t−1 .n, v, d + k +t − 1/ = Σi q̂i,t−1 .n, v, d + k + t − 1/. This amounts to rescaling the Z-marginals of the current
table as in step 4 of the algorithm (see below).
Now we shall write hij as hnvij to denote the dual solutions when the sample size is n andthe algorithm has run for v cycles. Similar updates will apply to the subscripts of β and γ.
Incorporating the above steps repeatedly over v complete cycles, and letting β1nvs = Σve=1 γ1nes ,β2nvr = Σve=1 γ2ner and β3nvt = Σve=1 γ3net , we obtainqij exp.hnvij /,qab exp.hnvab /q̂ij .n, v, d + k + k1 : : : kd / =.4:11/abwhere qij = q̂ij .n, 0, 0/ = .ni =n/mj =n = p̂ij , hnv = .hnvij / and k :::k

d
k
1 d
nr
mt
β1nvs .yi zsjs − cs / +β2nvr Iijr −β3nvt Jijt −+:hnvij =
n
n
s=1
r=1
t=1

.4:12/

Thus we propose the following algorithm for solving problem (4.1).
4.1. Algorithm
Step 1: calculate the constraint values c1 , : : : , cd and marginal totals ni =n and mj =n from theobserved table of .Y , Z/ values. Initialize as γn00 = 0 and qij = q̂ij .n, 0, 0/ = .ni =n/mj =n = p̂ij .
Let v = 1.

12

B. Bhattacharya and M. Al-talib

Step 2: begin with u = 1. At the vth cycle, uth step, let γ1nvu = γ Å solve equation (4.9). Useequation (4.10) to form q̂ij .n, v, u/. Do for u = 1, : : : , d.
Step 3: replace q̂ij .n, v, u/ by q̂ij .n, v, u/ni ={nq̂i+ .n, v, u/}, where q̂i+ .n, v, u/ = Σj q̂ij .n, v, u/,
1  i  k, d + 1  u  d + k.
Step 4: replace q̂ij .n, v, u/ by q̂ij .n, v, u/mj ={nq̂+j .n, v, u/}, where q̂+j .n, v, u/ = Σi q̂ij .n, v, u/,∀ j, d + k + 1  u  d + k + k1 : : : kd .
nv
nv
Step 5: Let q̂ij .n, v, d + k + k1 : : : kd / = q̂nv
ij . Stop if max{|Σij yi zsjs q̂ij − cs |, ∀ s, |q̂i+ − ni =n|,nv
∀i, |q̂+j − mj =n|, ∀ j} < , for some prespeciﬁed > 0; otherwise, replace v by v + 1, and go tostep 2.
We have used Fortran to implement the algorithm. This algorithm is for a ﬁxed sample sizen. Using results from Bhattacharya and Dykstra (1997), it follows that the above algorithm
converges as v → ∞. Convergence is in seconds for all the cases that we considered. Thus
hnv = .hnvij / converges to hn = .hnij / as v → ∞, where k :::k

d
k
1 d
nr
mt
hnij =
β1ns .yi zsjs − cs / +β2nr Iijr −β3nt Jijt −
+
,
n
n
s=1
r=1
t=1




d
nr
mt
=
β1ns .yi zsjs − cs / + β2ni 1 −+ β3nj 1 −nn
s=1for some β1ns , β2nr and β3nt . Comparing with equation (4.12), it must be that βnv1s → β1ns , βnv2r→ β2nr , βnv3t → β3nt , ∀ s, r, t as v → ∞. The asymptotic distribution properties of β1ns are studiedin Section 5.
5.

Asymptotic properties

For sample size n, at cycle v, deﬁne the vector β1nv = .β1nv1 , : : : , β1nvd /T . Next we show thatthe sequence of vectors {β1nv , v  1} which is generated by the algorithm converges to the true
parameter β1 = β in model (2.4) as v → ∞ and n → ∞.
Theorem 2. For a ﬁxed n, the algorithm in Section 4.1 yields the sequence of solutions{β1nv , v  1, n  1} at the vth cycle. Then β1nv converges to β1n as v → ∞, and β1n converges
to the true parameter β1 = β in model (2.4), with probability 1, as n → ∞.
Writing β1n = βn , in theorem 3 we consider the asymptotic normality properties of βn . Recall
that dQÅ =dP = qÅ .y, z/ is deﬁned in equation (2.3).

Theorem 3. Assume that QÅ exists and exp.h/dQÅ < ∞ for β in an open neighbourhoodof zero. Assume that ||yz||2 .dQÅ =dP/dQÅ < ∞. Then√Dn.βn − β/ → N.0, Σ−1 ΣÅ Σ−1 /,

where


ΣÅ =

.yz − c/.yz − c/T exp.h/dQÅ

,
dQÅ


Σ=

.5:1/
.yz − c/.yz − c/T exp.h/dP:
exp.h/dP

Minimum Relative Entropy Correlation Model

13

To estimate the asymptotic variance of βn , we replace the integrals in expression (5.1) withthe discrete sums, use midpoints of the discretized intervals as the values of yi and zij s, ﬁnd
inverses of the resulting matrices and use matrix multiplication.
6.

Simulation

6.1. Continuous case
Since the algorithm of Section 4.1 works for any distribution of Y and Zi s, we would like to compare its performance with some known procedure (e.g. maximum likelihood estimates (MLEs))
in cases when the correct model is known and when it is misspeciﬁed. We consider samples of sizesn = 50, 100, 150 for the random vector .Y , Z1 , Z2 / ∼ N3 .0, Σ/, where Σ = .2, 0:5, 0:5==0:5, 1:5,
0:5==0:5, 0:5, 2/. True β1 - and β2 -values are obtained by comparing the PDF of .Y , Z1 , Z2 / withequations (2.3) and (2.4) and computing Σ−1 ; here, true .β1 , β2 / = .0:1538, 0:1026/.
We need to categorize the observed values of Y and Z to construct their observed marginal
distributions. These are used to form the constraints in K. To determine the number of classes foreach variable, ﬁrst the optimal histogram bin width is calculated from Scott (1979), which asymptotically minimizes the integrated mean-squared error. He suggested the use of hn = 3:49σn−1=3 ,where σ is the standard deviation, when the sample is taken from√a normal
distribution. As the√√standard deviations of three variables are chosen to be close ( 2, 1:5, 2), we decided touse the same k for all three variables. We did some experimentation with different k for threevariables, but this did not produce any signiﬁcantly different results.
Using the values of y, z1 and z2 as the midpoints of those categories and frequency distributions, the observed marginal distributions of Y and .Z1 , Z2 / and the observed values of
E.YZ1 / and E.YZ2 / are found. These values are used to deﬁne the constraints in K. Multiplying
the observed marginal distributions of Y and .Z1 , Z2 /, a three-way contingency table is formed,
which is the starting point of the algorithm. The algorithm stops when all constraints are satisﬁed subject to a prespeciﬁed value. The algorithm converged in all cases. All simulations are
replicated 1000 times.
The results are in Table 1 using k = 8. A more detailed Table 1 is in the on-line supplement toshow the effect of the choice of k on the simulation results, where we varied k for n = 50, 100.
For instance, for k = 4, we considered intervals .−∞, − hn /, .hn , 0/, .0, hn / and .hn , ∞/, andsimilarly for other values of even k. The results are mostly invariant for k > 4. The results are
similar if the intervals are not centred near 0, or considering an odd number of intervals (thoseresults are not presented here).
The estimates that are produced by the algorithm are compared with the MLEs of β1 and β2 .
Comparing expression (2.3) with the expression of the normal PDF exponent .− 21 Σij xi σ ij xj /,−1
we ﬁnd that the MLEs of β1 and β2 are given by β̂ 1 = −σ̂ 12 and β̂ 2 = −σ̂ 13 , where Σ̂ = .σ̂ ij /is the inverse of the sample covariance matrix.
The misspeciﬁed model assumes that Y |Z ∼ exp.ZT β/ + N.0, 2/, Z ∼ N2 .0, ΣZ /, where ΣZ =
.1:5, 0:5==0:5, 2/ (the same as the covariance matrix of Z speciﬁed). This model is outside model(2.4) because the mean is exponential in β. But the procedure proposed is not affected as it does
not depend on Σ. See the on-line supplement for additional results from a wrongly speciﬁed
covariance matrix.
For each of β1 and β2 , we calculate the bias, standard error values and 95% coverage probability for β1 and β2 . Note that the MLE depends on the normal PDF, whereas our proceduredoes not. Simulations in Table 1 indicates that the procedure performs favourably with the MLE
when the correct model is used, and better than the MLE when the model is misspeciﬁed. The
performance of the procedure improves as the sample size increases.

14

B. Bhattacharya and M. Al-talib

Table 1.

Comparing the procedure with MLEs under correct and misspecified normal models

Methodnbias(β1n )se(β1n )bias(β2n )se(β2n )Procedure

50
100
150
50
100
150
50
100
150

−0:0065−0:0106−0:0099
0.01690.00970.0057−0:2908−0:2898−0:28900.09870.06440.05210.1116
0.07390.05850.10830.07030.0574−0:0091−0:0167−0:0137
0.00900.00310.0023−0:1896−0:1879−0:18990.08430.05530.0456
0.09630.06110.0503
0.0935
0.06120.0504MLE (correct)
MLE (misspeciﬁed)

CP(β1n )† CP(β2n )†
945
946
946
939
945
953
951
940
937

952
947
939
949
950
948
941
945
941

†95% coverage probability.

6.2. Discrete case
Suppose that Y takes values 0 and 1, with probabilities 0.4 and 0.6 respectively, and .Z1 , Z2 / =.0, 0/, .0, 1/, .1, 0/, .1, 1/ with probabilities 0.3, 0.25, 0.22 and 0.23 respectively. Assuming that
.β1 , β2 / = .1, 1/, an eight-cell multinomial probability vector is generated, where pi,j1 ,j2 ∝ ri,j1 ,j2 ×
exp.β1 yi z1j1 + β2 yi z2j2 + β3i + β4j1 j2 /, i, j1 , j2 = 1, 2, where the β3i - and β4j1 j2 -values are suchthat pi,j1 ,j2 satisﬁes the Y - and .Z1 , Z2 / marginal probabilities and ri,j1 ,j2 are obtained by multiplying the marginal probabilities of Y and .Z1 , Z2 /; ﬁnally, the true pi,j1 ,j2 -values are found bynumerical methods to be 0.1847, 0.0927, 0.0816, 0.0410, 0.1153, 0.1573, 0.1384 and 0.1890. The
constraints E.YZ1 / = c1 and E.YZ2 / = c2 simplify nicely to p221 + p222 = c1 and p212 + p222 =c2 respectively. We take random samples of size n from the eight-cell multinomial, for n =50, 100, 150. From the observed table we calculate the marginal distributions of Y and .Z1 , Z2 /,and also ﬁnd the observed values of c1 and c2 . In this case, exact dual solutions are available
not only for the marginal constraints but also for the correlation constraints; for any r = .rijk /,when solving minp:p221 +p222 =c1 I.p|r/, the dual problem is solved byγ̂ 1 = ln

c1 .r111 + r211 + r121 + r112 + r212 + r122 /,.1 − c1 /.r221 + r222 /and similarly for constraint p212 + p222 = c2 the dual problem is solved byγ̂ 2 = ln

c2 .r111 + r211 + r121 + r112 + r221 + r122 /:.1 − c2 /.r212 + r222 /As stated in Section 3.1, in the context of this simulation, we essentially have a logistic regression model. So we compared our procedure with the logistic regression model
P.Y = 1/ =exp.β0 + β1 yz1 + β2 yz2 /,1 + exp.β0 + β1 yz1 + β2 yz2 /where β0 , β1 and β2 are solved by using the Newton–Raphson (NR) method. As seen from
Table 2, the NR numbers are slightly lower than those from the procedure for n = 50; however,
it is the other way around for n = 100, 150. It is possible that this is due to the simple rescaling
nature of the algorithm that is proposed here; the NR method involves many matrix operations
including inversions, all of which contribute to the rounding errors, for example. The misspeciﬁed
model is taken to be a linear model pi,j1 ,j2 ∝ ri, j1 , j2 .β1 yi z1j1 + β2 yi z2j2 + β3i + β4j1 j2 /. Here theMinimum Relative Entropy Correlation Model
Table 2.
n
50
100
150

15

Comparing the procedure under correct and misspecified multinomial models
Methodbias(β1n )se(β1n )bias(β2n )se(β2n )Procedure
NR (correct)
Misspeciﬁed
Procedure
NR (correct)
Misspeciﬁed
Procedure
NR (correct)
Misspeciﬁed

0.11200.08203.17890.03980.03962.93010.03540.03982.91030.72850.67731.99090.44830.45001.29940.36910.37381.02410.07320.06532.41550.00750.02192.2793−0:0006
0.01542.20660.68420.6554
1.89460.45380.4553
1.20980.34630.34730.9828CP(β1n )† CP(β2n )†
954
946
696
953
953
411
951
952
195

950
946
783
951
950
551
955
957
391

†95% coverage probability.

probabilities are linear in β, so it is outside model (2.4). Here the procedure performs favourablywith the NR method when the model is correct, and it performs better when the model is
misspeciﬁed. See the on-line supplement for additional results when the model is misspeciﬁedwith β0 = 0 in the logistic regression model.
7.

Real data examples

We considered two examples with different correlation levels. Most of the covariates in the ﬁrst
and second example have respectively relatively low and high correlations with the response. We
ﬁt and compare different models in each case. The ﬁrst example has two continuous and one
discrete predictors, and it shows that the power transformation of the variables is useful. The
second example has four continuous predictors. In each case, the ﬁtted non-linear model seems
to be a better ﬁt to the data than the multiple linear regression model.
7.1. Trail making data
We consider a data set from Luo and Tsai (2012), which gives the scores of 334 patients on part
A of the trail making test with covariates education, age and diagnosis. The continuous response
variable Y is the test score, with range 0–150 s; if a patient cannot complete the test in 150 s, ascore of 150 is given. The continuous covariates are Z1 , years of education (3–21 years) and Z2 ,age (53–108 years), and the discrete covariate is Z3 , diagnosis (0 or 1). From the data, we ﬁndthat corr.Y , Z1 / = −0:40 corr.Y , Z2 / = 0:18 and corr.Y , Z3 / = 0:21.
To determine k, the number of classes, we ﬁnd, for each variable Z, the quantities sd, standard
deviation, and hn , optimal bin width (Scott, 1979), and then k is calculated as k = Z.334/ + −.Z.1/ − /=hn , where Z.334/ and Z.1/ are the maximum and minimum values of Z respectively
and > 0 is a very small number; of course, k is rounded to the nearest integer.
Using these classes, a four-way contingency table of sample proportions is formed fromthe observed data. Using this table and the midpoint of each class for each variable are used
to ﬁnd the observed values of E.YZ1 /, E.YZ2 / and E.YZ3 /, which deﬁne three dependence
constraints. In addition, observed Y -marginal totals produce 12 constraints and observed Zmarginal totals altogether produce 9 × 13 × 2 = 234 constraints, which make a total of 249
constraints. Multiplying the observed marginal distributions of Y and Z, a four-way contingency
table is formed, which is the starting point of the algorithm. The algorithm stops when all
constraints are satisﬁed subject to a prespeciﬁed value.

16

B. Bhattacharya and M. Al-talib

After the algorithm has converged, we ﬁnd model 1, β1Å = −4:76 × 10−3 , β2Å = 1:01 × 10−3 andÅ
β3 = 16:45 × 10−3 , with sd(β1Å / = 0:51 × 10−3 , sd(β2Å / = 0:09 × 10−3 and sd(β3Å / = 2:60 × 10−3 .
Thus all three coefﬁcients are statistically signiﬁcant. Although these values are on a par withthose in Luo and Tsai (2012), model 1 maintains the observed Y and Z marginals. Luo and Tsai(2012) found the interaction terms to be statistically insigniﬁcant, so we did not ﬁt them either.
If we keep Z2 and Z3 at a ﬁxed level, then the likelihood ratio that a patient’s score increases by
10 s will increase by a factor of exp{.0:00476/10} = 1:05, should the patient have 1 year less ofeducation, for all values of z1 .
The multiple linear regression line is obtained by using all data values and is given by ŷ =42:17 − 2:25z1 + 0:43z2 + 6:91z3 , with standard deviations of coefﬁcients 12.08, 0.32, 0.14 and2.41 respectively; thus all are signiﬁcant. Figs 1 and 2 show the plot of the ﬁtted model with nottransformed variables Z1 and Z2 (E.Y |Zi / versus Zi for i = 1, 2). The regression line is drawn by
ﬁxing the values of other predictors at their means. Also, 95% bootstrap bounds are shown ineach plot. Since our goal is to ﬁt a model based on dependence between Y and the covariates,
it would be of interest if the linear relationship between Y and some function of Zi s is any
stronger than that between Y and Zi s. Box and Tidwell (1962) suggested appropriate powertransformations of independent covariates that may be useful for our purpose. A plot of Z1
versus Z2 reveals that independence may be assumed (see the on-line supplement). Consideringαα.Y , Z1 1 , Z2 2 /, the Box–Tidwell method suggests the estimates α̂1 = −0:52 and α̂2 = 0. Hence,
we consider the transformations Z1−0:52 and log.Z2 /; note that corr.Y , Z1−0:52 / = −0:45 andcorr{Y , log.Z2 /} = 0:22 are slightly stronger than those without transformations. So we ﬁt the160
140

Trail Making A (s)

120
100
80
60
40
20
0

2

4

6

8

10

12

14

16

18

20

Years of Education
Fig. 1. Fitted trail making A (s) scores versus years of education using not transformed variables;
method proposed;, 95% bootstrap bounds;, multiple linear regression,

Minimum Relative Entropy Correlation Model

17

160
140

Trail Making A (s)

120
100
80
60
40
20
0

50

70

90

110

Age (years)
Fig. 2. Fitted trail making A (s) scores versus age (years) using not transformed variables;
proposed;, 95% bootstrap bounds;, multiple linear regression, methodmodel for Y by using the transformed covariates .Z1−0:52 , log.Z2 /, Z3 /. As before we choose kfor each variable in Table 4 (in the on-line supplement). A new four-way contingency table isformed by using the transformed variables, and similar steps are followed to those before. Then
we apply the algorithm for a total of 3 + 12 + 13 × 6 × 2 = 171 constraints, and ﬁnd model 2:β1ÅÅ = 0:25, β2ÅÅ = 0:07 and β3ÅÅ = 0:02, with sd(β1ÅÅ / = 45:32 × 10−3 , sd(β2ÅÅ / = 3:02 × 10−3 andsd(β3ÅÅ / = 3:18 × 10−3 . It is seen that, for model 2, all three coefﬁcients are statistically signiﬁcant
like in model 1. Next we also ﬁt a model using three interactions between the transformed
variables; however, the interactions turned out to be insigniﬁcant. So we stayed with the modelwith no interactions like in model 2.
Figs 3 and 4 show the plot of the ﬁtted model with transformed variables Z1 and Z2 . The same
regression line is drawn as before. Also, 95% bootstrap bounds are shown in each plot. Fromthese graphs, it is clear that model 2 attends to the non-linear nature of the data more than model
1 does. From equation (2.8), if we keep Z2 and Z3 at a ﬁxed level in model 2, then the likelihood
ratio that a patient’s score increases by 10 s will increase by a factor of exp[0:25{.z1 − 1/−0:5 −}10], should the patient have 1 year less education. Clearly this depends on the value of
z−0:5
1
z1 (as opposed to model 1), such as it equals 1.38, 1.06 and 1.01 when z1 = 3, 8, 21 respectively.
The residual plots that were obtained from the two models are given in Figs 9–12 in the on-line
supplement. From these plots, we observe that most residuals are around zero except a few
which have very large positive values. Those individuals performed much better than expected.
7.2. Prostate cancer data
Stamey et al. (1989) examined the correlations between the level of prostate-speciﬁc antigen andseveral clinical measures in 97 men who were about to receive a radical prostatectomy. The goal

18

B. Bhattacharya and M. Al-talib
160
140

Trail Making A (s)

120
100
80
60
40
20
0

2

4

6

8

10

12

14

16

18

20

Years of Education
Fig. 3.

Fitted trail making A (s) scores versus years of education using Box–Tidwell transformed variables;, method proposed;, 95% bootstrap bounds;, multiple linear regression

160
140

Trail Making A (s)

120
100
80
60
40
20
0

50

70

90

110

Age (years)
Fig. 4. Fitted trail making A (s) scores versus age (years) using Box–Tidwell transformed variables;
method proposed;, 95% bootstrap bounds;, multiple linear regression,

Minimum Relative Entropy Correlation Model

19(Hastie et al., 2009) is to predict Y , the logarithm of prostate-speciﬁc antigen level, lpsa, using
a number of measurements: Z1 , the logarithm of cancer volume, lcavol, Z2 , the logarithm of
prostate weight, lweight, Z3 , age, and Z4 , the logarithm of capsular penetration, lcp; all variables are continuous. We ﬁnd that corr.Y , Z1 / = 0:73, corr.Y , Z2 / = 0:43, corr.Y , Z3 / = 0:17 andcorr.Y , Z4 / = 0:55. Further, possible interactions between these covariates are also of interest,
such as, we ﬁnd that corr.Y , Z1 ÅZ2 / = 0:74, corr.Y , Z1 ÅZ3 / = 0:71 and corr.Y , Z1 ÅZ4 / = 0:36:
In this example, power functions of predictors did not signiﬁcantly increase correlations with Y .
The number of classes for the variables are determined as described earlier (see Table 4 in theon-line supplement). These classes are used to form a ﬁve-way contingency table of sample proportions from the observed data. Next, this table and the midpoints of each class for each variable
are used to ﬁnd the observed values of E.YZ1 /, E.YZ2 /, E.YZ3 / and E.YZ4 /, which deﬁne four
dependence constraints. In addition, the observed Y -marginal totals produce seven constraints
and observed Z-marginal totals altogether produce 6 × 4 × 7 × 7 = 1176 constraints, yielding a
total of 1187 constraints. The same steps are followed as in the previous example to ﬁnd theestimates of βi s. To decide on an appropriate model, ﬁrst we ﬁt the model with all four predictors.
(a) Model 1: predictors (lcavol, lweight, age, lcp), β̂ = .0:95, 1:40, −0:02, 0:37/ and sd.β̂/ =.0:08, 0:17, 0:01, 0:06/. Here ‘age’ is not a signiﬁcant predictor.
Next we consider the remaining three variables, with their interactions in model 2. The sameprocess as described above is followed.
(b) Model 2: predictors (lcavol, lweight, lcp, lcavol Å lweight, lcavol Å lcp, lweight Å lcp),β̂ = .1:02, 1:41, 0:40, −1:48, −0:57, −0:16/ and sd(β̂/ = .0:08, 0:18, 0:06, 0:87, 0:55, 0:38/.
Here none of the interaction terms is signiﬁcant. So we reﬁt the model without any of theinteraction terms with three predictors.
6
5
4

lpsa

3
2
1
0
-1
-1.5

-0.5

0.5

1.5

2.5

3.5

lcavolFig. 5. Fitted prostate lpsa-values versus lcavol from model 3;bootstrap bounds;, multiple linear regression, method proposed;

, 95%

20

B. Bhattacharya and M. Al-talib
6
5
4

lpsa

3
2
1
0
-1

2.3

2.8

3.3

3.8

4.3

4.8

lweight
Fig. 6. Fitted prostate lpsa-values versus lweight from model 3;bootstrap bounds;, multiple linear regression, method proposed;

, 95%

5.5
4.5

lpsa

3.5
2.5
1.5
0.5
-0.5
-1.5

-0.5

0.5

1.5

2.5

lcp
Fig. 7. Fitted prostate lpsa-values versus lcp from model 3;, method proposed;
strap bounds;, multiple linear regression (where lcp is not significant), 95% boot-(c) Model 3: predictors (lcavol, lweight, lcp), β̂= .0:91, 1:27, 0:38/ and sd.β̂/= .0:07, 0:06, 0:05/.
In Figs 5–7, we have plotted the conditional mean of lpsa given lcavol, lweight and lcp respectively, for model 3, and the 95% bootstrap bounds. The multiple linear regression curve is
given by lpsa = −0:73 + 0:58 lcavol + 0:67 lweight + 0:09 lcp, with standard deviations of thethree predictors 0:09, 0:18 and 0:07 respectively. Hence lcp is not a signiﬁcant predictor in theMinimum Relative Entropy Correlation Model21

linear regression model. In each ﬁgure, the regression line is drawn by ﬁxing the other predictors
at their mean values.
The residual plots that were obtained from model 3 are given in Figs 13–15 in the on-line
supplement. Since the residual plot versus lcp was deemed unsatisfactory, we have experimented
with values other than k = 4 for lcp and decided to go with k = 6 for lcp; however, this did notcause a big change in the values of β̂ i , and sd(β̂ i ), for any i. Thus all the results and ﬁguresthat are presented here use k = 6 for lcp. Scott (1979) mentioned that a value that is close to hissuggested value would be acceptable for data analysis.
8.

Discussion

This paper proposes a model that is characterized by the marginal distributions of the univariate
response, joint distribution of covariates and the correlations between the response and eachcovariate. A key feature of the method proposed is that no functional form on the conditional
distribution of the response given the covariates is assumed.
This paper uses the minimum relative entropy principle. An important point to note is thatusing this principle we ﬁnd a Q ∈ C which is closest to P (in KL distance measure), and safeguards
against using any other information in the process. Thus it may be that there are some very
speciﬁc deterministic relationships between Y and Z, but following this principle we might or
might not be able to detect it exactly because that speciﬁc relationship is not described in C.
Thus the optimality interpretation of theorem 1 is that the procedure proposed chooses thesolution from a pool (C) of all distributions that shares the same properties of correlation andmarginal distributions as those of the true (Y , Z) (but the exact relationship between Y and Z isunspeciﬁed) and is closest to the joint distribution of independent Y and Z.
To investigate the nature of the ﬁtted models for different possible relationships between Y andZ further, we consider three cases below. First we consider a straight line relationship between Y
and Z as in Y = 3 + 2Z, where Z ∼ unif.−1, 1/ and Y ∼ unif.3, 5/. Here corr.Y , Z/ = 1. We deﬁne
the set C as
C = {Q : E.YZ/ = c, Y ∼ P1 , Z ∼ P2 },where P1 = unif.3, 5/ and P2 = unif.−1, 1/, for some constant c. Let P = P1 P2 , and P1 and P2
be independent. Then P ∈ C. Although the exact I -projection of P onto C can be found as inequation (2.3), here we construct an estimate. On the basis of an observed sample of size 100
from the joint distribution of .Y , Z/, we form K as in equation (4.2); using the algorithm a ﬁtted
model is found, and this process is repeated 1000 times. Fig. 8(a) shows the actual model and theaverage of ﬁtted models from 1000 simulations. Here the averaged ﬁtted model matches exactly
the actual model.
Next we let Y and Z have a monotone relationship, but not a straight line. Let Y = exp.Z/={1 +exp.Z/}, where Z ∼unif.−5, 5/. Here corr.Y , Z/ ≈ 0:97, and corr.Y , Z3 / ≈ 0:8. Here Y ∼ fY .y/= 1={2y.1 − y/}, 0:007 < y < 0:993, and fY .y/ = 0 otherwise. Deﬁne C = {Q : E.YZ/ = c1 , E.YZ3 /= c2 , Y ∼ P1 , Z ∼ P2 }, where P1 corresponds to fY .y/ and P2 = unif.−5, 5/, for some constants
c1 and c2 . Let P = P1 P2 , and P1 and P2 be independent. Then the I -projection of P onto C is
the solution to the problem as in equation (2.3). To estimate, on the basis of an observed sampleof size 100, we form K. A ﬁtted model is found, and this process is repeated 1000 times. Fig. 8(b)shows the actual model and the average of ﬁtted models from 1000 simulations. Our procedure
cannot locate the speciﬁc relationship Y = exp.Z/={1 + exp.Z/} exactly, but it can ﬁnd a model
which satisﬁes the properties described in C, and is closest to P. The ﬁtted model would becloser to the actual model if more correlation constraints are considered, e.g. E.YZ5 / = c3 , in C.

22

B. Bhattacharya and M. Al-talib
1

6

0.9
0.8

5

0.7

4

0.6

y 0.5

y 3

0.4

2

0.3
0.2

1
0

0.1

1

2

3

0

-5

-3.75

-1.25

1.25

3.75

5

z
(b)

z
(a)
1
0.9
0.8
0.7
0.6

y 0.5
0.4
0.3
0.2
0.1
0

-1

-0.8 -0.6 -0.4 -0.2

0

0.2

0.4

0.6

0.8

1

z
(c)

Fig. 8. (a) True and fitted models coincide, (b) monotone relationship (, true model;
model) and (c) non-monotone relationship (, true model;, fitted model), fitted

Next suppose that the relationship between Y and Z is not monotone. Let Y = |Z| whereZ ∼unif.−1, 1/. Here corr.Y , Z/ = 0. Let
C = {Q : E.YZ/ = 0, Y ∼ P1 , Z ∼ P2 },where P1 = unif.0, 1/ and P2 = unif.−1, 1/. Let P = P1 P2 and P1 and P2 be independent. ThenP ∈ C, and the I -projection of P onto C is itself, or β = 0 in the expression (2.4) of the joint
distribution of the solution. Thus we cannot locate the true model.
However, note that corr.Y , Z2 / ≈ 0:97. Incorporating this information, deﬁneC = {Q : E.YZ/ = 0, E.YZ2 / = c, Y ∼ P1 , Z ∼ P2 },for some constant c. Then P ∈ C (C ⊂ C). Again, the I -projection of P onto C is given byequation (2.3). To estimate, on the basis of an observed sample of size 100, we form K anda ﬁtted model is found. This is repeated 1000 times. Fig. 8(c) shows the actual model and theMinimum Relative Entropy Correlation Model

23

average of ﬁtted models from 1000 simulations. Here the averaged ﬁtted model is close to theactual model.
It may be concluded that correlations alone may not be sufﬁcient to characterize the underlying dependence structure completely. However, the last example shows that including theadditional constraint corr(Y , Z2 ) in C is helpful in ﬁnding a suitable ﬁtted model. Thus werecommend power transformations of Y and/or Z in constructing C. See also Section 7.
The sensitivity of the choice of transformations can be seen in the trail making data example
where we presented results with and without transformations. We have proposed the use of the
Box and Tidwell (1962) method to ﬁnd a suitable transformation, which requires independent
covariates. For real data with dependent covariates, orthogonality of covariate vectors may be−1=2
obtained by using Sz Z as covariates instead of Z, where Sz = cov(Z).
The characterization of known models like those described in Section 3 through dependence
using correlations between the response and each covariate appears to be new. This observation
has led to the algorithm that is given in the paper. The practical advantage is that the algorithm
proposed always converges, is simple to use (e.g. no matrix inversion as in the NR method and isapplicable for any response and covariates. Robustness of the method to model misclassiﬁcation
is demonstrated in simulations for both continuous and discrete cases. The reason for this beneﬁt
is that it does not depend on any likelihood, as opposed to maximum likelihood estimators.
Acknowledgements
The authors thank the Joint Editor, the Associate Editor and the referees for their detailed
reviews and encouraging comments. They helped to improve the manuscript greatly.
References
Agresti, A. (2013) Categorical Data Analysis. New York: Wiley.
Bhattacharya, B. (2006) An iterative procedure for general probability measures to obtain I-projection onto
intersection of convex sets. Ann. Statist., 34, 878–902.
Bhattacharya, B. and Dykstra, R. (1995) A general duality approach to I-projections. J. Statist. Planng Inf., 47,
203–216.
Bhattacharya, B. and Dykstra, R. (1997) A Fenchel duality aspect of iterative I-projection procedures. Ann. Inst.
Statist. Math., 49, 435–446.
Box, G. E. P. and Tidwell, P. W. (1962) Transformation of the independent variables. Technometrics, 4, 531–542.
Chan, K. C. G. (2013) Nuisance parameter elimination for proportional likelihood ratio models with nonignorable
missingness and random truncation. Biometrika, 100, 269–276.
Csiszár, I. (1975) I-divergence geometry of probability distributions and minimization problems. Ann. Probab., 3,
146–158.
Fokianos, K. and Kaimi, I. (2006) On the effect of misspecifying the density ratio model. Ann. Inst. Statist. Math.,
58, 475–497.
Gilbert, P. B., Lele, S. R. and Vardi, Y. (1999) Maximum likelihood estimation in semiparametric selection bias
models with application to AIDS vaccine trials. Biometrika, 86, 27–43.
Good, I. J. (1976) On the application of symmetric Dirichlet distributions and their mixtures to contingency tables.
Ann. Statist., 4, 1159–1189.
Guererro, V. M. and Johnson, R. A. (1982) Use of the Box-Cox transformation with binary response models.
Biometrika, 69, 309–314.
Hastie, T., Tibshirani, R. and Friedman, J. (2009) The Elements of Statistical Learning: Data Mining, Inference,
and Prediction, 2nd edn. New York: Springer.
Huang, A. and Rathouz, P. J. (2012) Proportional likelihood ratio models for mean regression. Biometrika, 99,
223–229.
Ichimura, H. (1993) Semiparametric least squares (SLS) and weighted SLS estimation of single-index models. J.
Econmetr., 58, 71–120.
Jaynes, E. T. (1957) Information theory and statistical mechanics. Phys. Rev., 106, 620–630.
Joe, H. (1989) Relative entropy measures of multivariate dependence. J. Am. Statist. Ass., 84, 157–164.
Kay, R. and Little, S. (1987) Transformations of the explanatory variables in the logistic regression model forbinary data. Biometrika, 74, 495–501.

24

B. Bhattacharya and M. Al-talib

Luo, X. and Tsai, W. Y. (2012) A proportional likelihood ratio model. Biometrika, 99, 211–222.
Qin, J. (1998) Inferences for case-control data and semiparametric two-sample density ratio models. Biometrika,
85, 619–630.
Rathouz, P. J. and Gao, L. P. (2009) Generalized linear models with unspeciﬁed reference distribution. Biostatistics,
10, 205–218.
Reshef, D. N., Reshef, Y. A., Finucane, H. G., Grossman, S. R., McVean, G., Turnbaugh, P. J., Lander, E.
S., Mitzenmacher, M. and Sabeti, P. C. (2011) Detecting novel associations in large data sets. Science, 334,
1518–1522.
Rockafellar, R. T. (1970) Convex Analysis. Princeton: Princeton University Press.
Scott, W. D. (1979) On optimal and data-based histograms. Biometrika, 66, 605–610.
Stamey, T. A., Kabalin, J. N., McNeal, J. E., Johnstone, I. M., Freiha, F., Redwine, E. A. and Yang, N. (1989)
Prostate speciﬁc antigen in the diagnosis and treatment of adenocarcinoma of the prostate: II, Radical prostatectomy treated patients. J. Urol., 141, 1076–1083.
Szekely, G. J. and Rizzo, M. L. (2009) Brownian distance covariance. Ann. Appl. Statist., 3, 1236–1265.
Szekely, G. J., Rizzo, M. L. and Bakirov, N. K. (2007) Measuring and testing dependence by correlation of
distances. Ann. Statist., 35, 2769–2794.

Supporting information
Additional ‘supporting information’ may be found in the on-line version of this article:‘Supplement to: A minimum relative entropy based correlation model between the response and covariates’.