On the role of the overall effect in exponential families

TL;DR

This paper examines how adding or removing the overall effect in exponential families affects their properties, geometry, and computational aspects, with implications for statistical modeling and biological data analysis.

Contribution

It characterizes the impact of the overall effect on exponential family properties, geometry, and algorithms, linking algebraic geometry concepts to statistical modeling.

Findings

Adding the overall effect creates the smallest regular exponential family containing the curved one.

Removing the overall effect simplifies the family but can lead to different estimation properties.

Including the overall effect can produce estimates outside the intended model in biological applications.

Abstract

Exponential families of discrete probability distributions when the normalizing constant (or overall effect) is added or removed are compared in this paper. The latter setup, in which the exponential family is curved, is particularly relevant when the sample space is an incomplete Cartesian product or when it is very large, so that the computational burden is significant. The lack or presence of the overall effect has a fundamental impact on the properties of the exponential family. When the overall effect is added, the family becomes the smallest regular exponential family containing the curved one. The procedure is related to the homogenization of an inhomogeneous variety discussed in algebraic geometry, of which a statistical interpretation is given as an augmentation of the sample space. The changes in the kernel basis representation when the overall effect is included or removed are derived. The geometry of maximum likelihood estimates, also allowing zero observed frequencies, is described with and without the overall effect, and various algorithms are compared. The importance of the results is illustrated by an example from cell biology, showing that routinely including the overall effect leads to estimates which are not in the model intended by the researchers.

Equations

$$RM(\mathbf{A})=\{\mbox{\boldmath$p$}\in\mathcal{P}:\,\,\log\mbox{\boldmath$p$}=\mathbf{A}^{\prime}\mbox{\boldmath$\theta$}\},$$

$$RM(\mathbf{A})=\{\mbox{\boldmath$p$}\in\mathcal{P}:\,\,{\mathbf{D}}\log\mbox{\boldmath$p$}=\mbox{\boldmath$0$}\}.$$

$${\mathbf{p}}^{{\mathbf{d}}_1^{+}}/{\mathbf{p}}^{{\mathbf{d}}_1^{-}}=1,{\mathbf{p}}^{{\mathbf{d}}_2^{+}}/{\mathbf{p}}^{{\mathbf{d}}_2^{-}}=1,\cdots {\mathbf{p}}^{{\mathbf{d}}_K^{+}}/{\mathbf{p}}^{{\mathbf{d}}_K^{-}}=1,$$

$${\mathbf{p}}^{{\mathbf{d}}_1^{+}}-{\mathbf{p}}^{{\mathbf{d}}_1^{-}}=0,{\mathbf{p}}^{{\mathbf{d}}_2^{+}}-{\mathbf{p}}^{{\mathbf{d}}_2^{-}}=0,\cdots {\mathbf{p}}^{{\mathbf{d}}_K^{+}}-{\mathbf{p}}^{{\mathbf{d}}_K^{-}}=0,$$

$${\mathcal{X}}_{\mathbf{A}}=\left\{\mathbf{p}\in {\mathbb{R}}_{\ge 0}^{|\mathcal{I}|}:{\mathbf{p}}^{{\mathbf{d}}^{+}}={\mathbf{p}}^{{\mathbf{d}}^{-}},\forall \mathbf{d}\in Ker(\mathbf{A})\right\}.$$

$$RM(\mathbf{A})={\mathcal{X}}_{\mathbf{A}}\cap int({\Delta }_{I-1}),$$

$$\text{vbox}\text{hrule}height=0.5ptRM(\mathbf{A})=\mathbf{p}\in {\mathcal{X}}_{\mathbf{A}}\cap {\Delta }_{I-1}.$$

$$\mbox{\boldmath$d$}_{1}^{\prime}\mbox{\boldmath$1$}\neq 0,\quad\mbox{\boldmath$d$}_{2}^{\prime}\mbox{\boldmath$1$}=0,\quad\dots,\,\,\mbox{\boldmath$d$}_{K}^{\prime}\mbox{\boldmath$1$}=0.$$

$$\bar{\mathbf{A}}=\left(\begin{array}[]{c}\boldsymbol{1}^{\prime}\\ \mathbf{A}\end{array}\right).$$

$$\mbox{\boldmath$d$}_{1}^{\prime}\log\mbox{\boldmath$p$}=0,\quad\mbox{\boldmath$d$}_{2}^{\prime}\log\mbox{\boldmath$p$}=0,\dots,\quad\mbox{\boldmath$d$}_{K}^{\prime}\log\mbox{\boldmath$p$}=0,\quad\mbox{ for any }\,\,\boldsymbol{p}\in RM({\mathbf{A}}).$$

$$\bar{\mathbf{D}}=\left(\begin{array}[]{c}\boldsymbol{d}_{2}^{\prime}\\ \vdots\\ \boldsymbol{d}_{K}^{\prime}\end{array}\right).$$

$$\bar{\mathbf{D}}\bar{\mathbf{A}}^{\prime}=\left(\begin{array}[]{c}\boldsymbol{d}_{2}^{\prime}\\ \vdots\\ \boldsymbol{d}_{K}^{\prime}\end{array}\right)\left(\begin{array}[]{cc}\boldsymbol{1}&{\mathbf{A}}^{\prime}\end{array}\right)=\left(\begin{array}[]{cc}\boldsymbol{d}_{2}^{\prime}\boldsymbol{1}&\boldsymbol{d}_{2}^{\prime}\mathbf{A}^{\prime}\\ \vdots&\vdots\\ \boldsymbol{d}_{K}^{\prime}\boldsymbol{1}&\boldsymbol{d}_{K}^{\prime}\mathbf{A}^{\prime}\end{array}\right)=\mathbf{O},$$

$$\mbox{\boldmath$d$}_{2}^{\prime}\log\mbox{\boldmath$p$}=0,\dots,\quad\mbox{\boldmath$d$}_{K}^{\prime}\log\mbox{\boldmath$p$}=0,\quad\mbox{ for any }\,\,\boldsymbol{p}\in RM(\bar{\mathbf{A}}).$$

$$\log p_i={\theta }_0+a_{1i}{\theta }_1+\cdots +a_{Ji}{\theta }_J,$$

$$RM(\mathbf{A})=\{\mathbf{p}\in RM(\bar{\mathbf{A}}):{\theta }_0(\mathbf{p})=0\}.$$

$$\mathbf{D}\log\boldsymbol{p}=\left(\begin{array}[]{c}\mbox{\boldmath$d$}_{1}\\[3.0pt] \bar{\mathbf{D}}\end{array}\right)\log\boldsymbol{p}=\left(\begin{array}[]{c}\theta_{0}(\boldsymbol{p})\mbox{\boldmath$d$}_{1}^{\prime}\mbox{\boldmath$1$}\\[3.0pt] \bar{\mathbf{D}}\mathbf{A}^{\prime}\boldsymbol{\theta}\end{array}\right)=\left(\begin{array}[]{c}-\theta_{0}(\boldsymbol{p})\\[3.0pt] \mbox{\boldmath$0$}\end{array}\right).$$

$$\Theta=\{\mbox{\boldmath$\theta$}^{\prime}=(\theta_{1},\dots,\theta_{J})^{\prime}\in\mathbb{R}^{J}:\,\,\mbox{\boldmath$1$}^{\prime}\exp\{\mathbf{A}^{\prime}\mbox{\boldmath$\theta$}\}=1\}.$$

$$\Theta_{1}=\{(\theta_{0},\mbox{\boldmath$\theta$}^{\prime})=(\theta_{0},\theta_{1},\dots,\theta_{J})^{\prime}\in\mathbb{R}^{J+1}:\,\,\theta_{0}=-\log(\mbox{\boldmath$1$}^{\prime}\exp\{\mathbf{A}^{\prime}\mbox{\boldmath$\theta$}\})\}.$$

$$\mathbf{A}=\left(\begin{array}[]{rrrr}1&1&1&0\\ 0&0&1&1\\ \end{array}\right),\qquad\bar{\mathbf{A}}=\left(\begin{array}[]{rrrr}1&1&1&1\\ 1&1&1&0\\ 0&0&1&1\\ \end{array}\right)$$

$$\left\{\begin{array}[]{l}p_{1}=\alpha_{1},\\ p_{2}=\alpha_{1},\\ p_{3}=\alpha_{1}\alpha_{2},\\ p_{4}=\alpha_{2},\end{array}\right.\qquad\left\{\begin{array}[]{l}p_{1}=\beta_{0}\beta_{1},\\ p_{2}=\beta_{0}\beta_{1},\\ p_{3}=\beta_{0}\beta_{1}\beta_{2},\\ p_{4}=\beta_{0}\beta_{2},\end{array}\right.$$

$$\left\{\begin{array}[]{l}\mbox{\boldmath$d$}_{1}^{\prime}\log\mbox{\boldmath$p$}=0,\\ \mbox{\boldmath$d$}_{2}^{\prime}\log\mbox{\boldmath$p$}=0,\end{array}\right.\qquad\left\{\begin{array}[]{l}\mbox{\boldmath$d$}_{2}^{\prime}\log\mbox{\boldmath$p$}=0.\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}p_{3}/(p_{1}p_{4})=1,\\ p_{1}/p_{2}=1,\end{array}\right.\qquad\left\{\begin{array}[]{l}p_{1}/p_{2}=1.\\ \end{array}\right.$$

$$\bar{\mathbf{A}}_{1}=\left(\begin{array}[]{cc}\boldsymbol{1}_{(I-I_{0})}^{\prime}&\boldsymbol{1}_{I_{0}}^{\prime}\\ \mathbf{A}_{1}&\mathbf{O}_{(J-1)\times I_{0}}\end{array}\right).$$

$${\mathbf{d}}_i=({\mathbf{u}}_i^{\prime },{\mathbf{v}}_i^{\prime }{)}^{\prime },\text{mbox}fori=1,\dots ,I-J.$$

$$\boldsymbol{0}=\bar{\mathbf{A}}_{1}\boldsymbol{d}_{i}=\left(\begin{array}[]{cc}\boldsymbol{1}_{(I-I_{0})}^{\prime}&\boldsymbol{1}_{I_{0}}^{\prime}\\ \mathbf{A}_{1}&\mathbf{O}\end{array}\right)\left(\begin{array}[]{c}\boldsymbol{u}_{i}\\ \boldsymbol{v}_{i}\end{array}\right),\quad\mbox{ for }\quad i=1,\dots,I-J,$$

$$1_{(I-I_0)}^{\prime }{\mathbf{u}}_i+1_{I_0}^{\prime }{\mathbf{v}}_i=0,{\mathbf{A}}_1{\mathbf{u}}_i=0,\text{mbox}fori=1,\dots ,I-J.$$

$${\mathbf{d}}_1$$

$${\mathbf{d}}_2$$

$$\cdots$$

$${\mathbf{d}}_{I_0}$$

$${\mathbf{d}}_{I_0+1}$$

$${\mathbf{d}}_{I_0+2}$$

$$\cdots$$

Full text

• Keywords:

algebraic variety, contingency table, independence, log-linear model, maximum likelihood estimation, overall effect, relational model

1 Introduction

This paper deals with exponential families of probability distributions over discrete sample spaces. When defining such families, usually, a normalizing constant, which of course, is constant over the sample space but not over the family, is included. The presence of the normalizing constant implies that the parameter space may be an open set, which, in turn, is necessary for asymptotic normality of estimates and for the applicability of standard testing procedures. The normalizing constant, from an applied perspective, may be interpreted as a baseline or common effect, present everywhere on the sample space and is, therefore, also called the overall effect. The focus of the present work is to better understand the implications of having or not having an overall effect in such families, in particular how adding or removing the overall effect affects the properties of discrete exponential families.

Motivated by a number of important applications, Klimova \BOthers. (\APACyear2012); Klimova \BBA Rudas (\APACyear2012, \APACyear2016) developed the theory of relational models, which generalize discrete exponential families, also called log-linear models, to situations when the sample space is not necessarily a full Cartesian product, the statistics defining the exponential family are not necessarily indicators of cylinder sets, and the overall effect is not necessarily present. Exponential families without the overall effect are particularly relevant, sometimes necessary, when the sample space is a proper subset of a Cartesian product. Several real examples, when certain combinations of the characteristics were either not possible logically or were left out from the design of the experiment were discussed in Klimova \BOthers. (\APACyear2012). A real problem of this structure from cell biology is analyzed in this paper, too. When the overall effect is not present, the standard normalization procedure to obtain probability distributions cannot be applied, because the family is curved Klimova \BOthers. (\APACyear2012). When, in spite of this, the standard normalization procedure is applied, as was done in this analysis, the resulting estimates do not possess the fundamental model properties.

The standardization of the estimates in exponential families is also an issue, when the size of the problem is very large and the computational burden is significant. Some Neural probabilistic language models are relational models. Due to the high-dimensional sample space, the evaluation of the partition function, which is needed for normalization, may be intractable. Some of the methods of parameter estimation under such models are based on the removal of the partition function, that is, the removal of the overall effect from the model and performing model training using the models without the overall effect. Approximations of estimates with and without the overall effect were studied, for example, by Mnih \BBA Teh (\APACyear2012) and Andreas \BBA Klein (\APACyear2015), among others. A different approach to avoiding global normalization (i.e., having an overall effect) is described in Koller \BBA Friedman (\APACyear2009). However, the implications of the removal of the overall effect are not discussed in the existing literature.

Another area where removing or including the overall effect is relevant, is context specific independence models, see, e.g., Høsgaard (\APACyear2004) and Nyman, Pensar, Koski\BCBL \BBA Corander (\APACyear2016). When the sample space is an incomplete Cartesian product, removing the overall effect, as described in this paper, specifies different variants of conditional independence in the parts of the sample space, depending on whether or not the part is or is not affected affected by the missing cells.

While including the overall in the definition of the statistical model to be investigated is seen by many researchers as “natural’ or “harmless”, we show in this paper that adding or removing the overall effect may dramatically change the characteristics of the exponential family, up to the point of altering the fundamental model property intended by the researcher.

The main results of the paper include showing that allowing the overall effect expands the curved exponential family to the smallest regular exponential family which contains it. When the overall effect is removed, the sample space may have to be reduced (if there were cells which contained the overall effect only), and the changes in the structure of the generalized odds ratios defining the model are described in both cases. In the language of algebraic geometry, the procedure of removing the overall effect is identical to the dehomogenization of the variety defining the model (Cox, Little\BCBL \BBA O’Shea, \APACyear2007). An important area of applications of the results presented here is the case when several binary features are observed, but the combination that no feature is present is either is impossible logically or is possible but was left out from the study design. The converse of dehomogenization, that is homogenizing a variety, involves including a new variable, and it is shown that in some cases this can be identified, from a statistical perspective, with augmenting the sample space by a cell which is characterized by no feature being present. For example, the Aitchison – Silvey independence (Aitchison \BBA Silvey, \APACyear1960; Klimova \BBA Rudas, \APACyear2015) is homogenized, through the augmentation of the sample space, into the standard independence model.

The paper is organized as follows. Section 2 gives a canonical definition of relational models using homogeneous, and if there is no overall effect included, one inhomogeneous generalized odds ratios, called dual representation and shows that including the overall effect is identical to omitting the inhomogeneous generalized odds ratio from it.

Section 3 contains the result that including the overall effect expands the curved exponential family into the smallest regular one containing it. For the case of the removal of the overall effect, the dual representation of the model is given, and the relevance of certain results in algebraic geometry to the statistical problem is discussed. In particular, the homogenization of a variety through including a new variable is identified with augmenting the sample space with a cell where no feature is present, when this is meaningful. It is proved that the homogenization of the Aitchison – Silvey (in the sequel, AS) independence model, which is defined on sample spaces where all combinations of features, except for the “no feature present” combination, are possible, is the usual model of mutual independence on the full Cartesian product obtained after augmenting the sample space with the missing cell. The relationship of these results with context specific independence is also described.

Section 4 compares the maximum likelihood (ML) estimates in geometrical terms for relational models with and without the overall effect and based on the insight obtained, a modification of the algorithm proposed in Klimova \BBA Rudas (\APACyear2015) is given. It is illustrated, that the ML estimates under two models which differ only in the lack or presence of the overall term, may be very different, up to the point of the existence or no existence of positive ML estimates, when the data contain observed zeros. However, when the MLE exists in the model containing the overall effect, it also does in the model obtained after the removal of the overall effect.

Finally, Section 5 discusses an example of applications of relational models in cell biology. The equal loss of potential model in hematopoiesis (Perié \BOthers., \APACyear2014) is a relational model without the overall effect. The published analysis of this model added the overall effect to it, to simplify calculations, and with this changed the properties of the model so that the published estimates do not fulfill the fundamental model property.

2 A canonical form of relational models

Let $Y_{1},\dots,Y_{K}$ be random variables taking values in finite sets $\mathcal{Y}_{1},\dots,\mathcal{Y}_{K}$, respectively. Let the sample space $\mathcal{I}$ be a non-empty, proper or improper, subset of $\mathcal{Y}_{1}\times\dots\times\mathcal{Y}_{K}$, written as a sequence of length $I=|\mathcal{I}|$ in the lexicographic order. Assume that the population distribution is parameterized by cell probabilities $\boldsymbol{p}=(p_{1},\dots,p_{I})$, where $p_{i}\in(0,1)$ and $\sum_{i=1}^{I}p_{i}=1$, and denote by $\mathcal{P}$ the set of all strictly positive distributions on $\mathcal{I}$. For simplicity of exposition, a distribution in $\mathcal{P}$ will be identified with its parameter, $\boldsymbol{p}$, and \mathcal{P}=\{\mbox{\boldmathp}>\mbox{\boldmath0}:\,\,\mbox{\boldmath1}^{\prime}\mbox{\boldmathp}=1\}.

Let $\mathbf{A}$ be a $J\times I$ matrix of full row rank with 0–1 elements and no zero columns. A relational model for probabilities $RM(\mathbf{A})$ generated by $\mathbf{A}$ is the subset of $\mathcal{P}$ that satisfies:

[TABLE]

where $\theta$ = $(\theta_{1},\dots,\theta_{J})^{\prime}$ is the vector of log-linear parameters of the model. A dual representation of a relational model can be obtained using a matrix, $\mathbf{D}$, whose rows form a basis of $Ker(\mathbf{A})$, and thus, $\mathbf{D}\mathbf{A}^{\prime}$ = $\mathbf{O}$:

[TABLE]

The number of the degrees of freedom $K$ of the model is equal to $dim(Ker(\mathbf{A}))$. In the sequel, \mbox{\boldmathd}_{1}^{\prime},\mbox{\boldmathd}_{2}^{\prime},\dots,\mbox{\boldmathd}_{K}^{\prime} denote the rows of $\mathbf{D}$. The dual representation can also be expressed in terms of the generalized odds ratios:

[TABLE]

or in terms of the cross-product differences:

[TABLE]

where ${\boldsymbol{d^{+}}}$ and ${\boldsymbol{d^{-}}}$ stand for, respectively, the positive and negative parts of a vector $\boldsymbol{d}$. The following dual representation is invariant of the choice of the kernel basis.

Let $\mathcal{X}_{\mathbf{A}}$ denote the polynomial variety associated with $\mathbf{A}$ (Sturmfels, \APACyear1996):

[TABLE]

The relational model generated by $\mathbf{A}$ is the following set of distributions:

[TABLE]

where $int(\Delta_{I-1})$ is the interior of the $(I-1)$-dimensional simplex.

Notice that the variety $\mathcal{X}_{\mathbf{A}}$ includes elements $\boldsymbol{p}$ with zero components as well and can be used to extend the definition of the model to allow zero probabilities. The extended relational model, \hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-2.5ptRM\kern-1.00006pt}}}(\mathbf{A}), is the intersection of the variety $\mathcal{X}_{\mathbf{A}}$ with the probability simplex:

[TABLE]

See Klimova \BBA Rudas (\APACyear2016) for more detail on the extended relational models.

Let \mbox{\boldmath1}^{\prime} = $(1,\dots,1)$ be the row of $1$’s of length $I$. If $\boldsymbol{1}^{\prime}$ does not belong to the space spanned by the rows of $\mathbf{A}$, the relational model $RM(\mathbf{A})$ is said to be a model without the overall effect. Such models are specified using homogeneous and at least one non-homogeneous generalized odds ratios, and the corresponding variety $\mathcal{X}_{\mathbf{A}}$ is non-homogeneous (Klimova \BOthers., \APACyear2012).

Proposition 1**.**

Let $RM(\mathbf{A})$ be a model without the overall effect. There exists a kernel basis matrix $\mathbf{D}$ whose rows satisfy:

[TABLE]

Proof.

A relational model does not contain the overall effect if and only if it can be written using non-homogeneous (and possibly homogeneous) generalized odds ratios (Klimova \BOthers., \APACyear2012). Therefore, $\mathbf{D}$ has at least one row, say \mbox{\boldmathd}_{1}^{\prime}, that is not orthogonal to $1$: C_{1}=\mbox{\boldmathd}_{1}^{\prime}\mbox{\boldmath1}\neq 0.

Suppose there exists another row, say $\boldsymbol{d}_{2}^{\prime}$, that is not orthogonal to $\boldsymbol{1}$ and thus $C_{2}=\boldsymbol{d}_{2}^{\prime}\boldsymbol{1}\neq 0$. The vectors $\boldsymbol{d}_{1}$ and $\boldsymbol{d}_{2}$ are linearly independent, so are the vectors $\boldsymbol{d}_{1}$ and $C_{2}\boldsymbol{d}_{1}-C_{1}\boldsymbol{d}_{2}$. Substitute the row $\boldsymbol{d}_{2}^{\prime}$ with the row $C_{2}\boldsymbol{d}_{1}^{\prime}-C_{1}\boldsymbol{d}_{2}^{\prime}$. And so on. ∎

It is assumed in the sequel that $\boldsymbol{1}^{\prime}$ is not in the row space of $\mathbf{A}$. Notice that, because $\mathbf{A}$ is 0–1 matrix without zero columns, this is only possible when $2\leq J=rank(\mathbf{A})<I-1$. Throughout the entire paper, the kernel basis matrix $\mathbf{D}$ is assumed to satisfy (8), and, without loss of generality, \mbox{\boldmathd}_{1}^{\prime}\mbox{\boldmath1}=-1.

Some consequences of adding the overall effect to a relational model will be investigated by comparing the properties of the relational model generated by $\mathbf{A}$ and the model generated by the matrix $\bar{\mathbf{A}}$ obtained by augmenting the model matrix $\mathbf{A}$ with the row $\boldsymbol{1}^{\prime}$:

[TABLE]

Let $RM(\bar{\mathbf{A}})$ be the relational model generated by $\bar{\mathbf{A}}$. Because $\boldsymbol{1}^{\prime}$ is a row of $\bar{\mathbf{A}}$, the corresponding polynomial variety $\mathcal{X}_{\bar{\mathbf{A}}}$ is homogeneous (cf. Sturmfels, \APACyear1996).

Theorem 1**.**

The dual representation of $RM(\bar{\mathbf{A}})$ can be obtained from the dual representation of $RM({\mathbf{A}})$ by removing the constraint specified by a non-homogeneous odds ratio from the latter.

Proof.

Write the dual representation of $RM({\mathbf{A}})$ in terms of the generalized log odds ratios:

[TABLE]

By the previous assumption, \mbox{\boldmathd}_{1}^{\prime}\mbox{\boldmath1}=-1, and thus, the constraint \mbox{\boldmathd}_{1}^{\prime}\log\mbox{\boldmathp}=0 is specified by a non-homogeneous odds ratio. Define $\bar{\mathbf{D}}$ as:

[TABLE]

Because, from (8), \mbox{\boldmathd}_{2},\dots,d_{K}\in Ker(\mathbf{A}),

[TABLE]

and thus, \mbox{\boldmathd}_{2},\dots,d_{K}\in Ker(\bar{\mathbf{A}}). Finally, as $rank(\bar{\mathbf{D}})=K-1$, \mbox{\boldmathd}_{2},\dots,d_{K} is a basis of $Ker(\bar{\mathbf{A}})$, and therefore,

[TABLE]

∎

3 The influence of the overall effect on the model structure

The consequences of adding or removing the overall effect will be studied separately. The changes in the model structure after the overall effect is added are considered first.

Let $RM(\mathbf{A})$ be a relational model without the overall effect and $RM(\bar{\mathbf{A}})$ be the corresponding augmented model. Let $\mathbf{A}=(a_{ji})$ for $j=1,\dots,J$, $i=1,\dots,I$. For any $\boldsymbol{p}\in RM(\bar{\mathbf{A}})$:

[TABLE]

where $\theta_{j}=\theta_{j}(\boldsymbol{p})$, $j=0,1,\dots,J$, are the log-linear parameters of $\boldsymbol{p}$. In particular, $\theta_{0}(\boldsymbol{p})$ is the overall effect of $\boldsymbol{p}$.

Theorem 2**.**

The augmented model, $RM(\bar{\mathbf{A}})$, is the minimal regular exponential family which contains $RM(\mathbf{A})$, and

[TABLE]

Proof.

The second claim is proved first. Denote $\mathcal{M}_{0}=\{\boldsymbol{p}\in RM(\bar{\mathbf{A}}):\,\,\theta_{0}(\boldsymbol{p})=0\}$. Let $\mathbf{D}$ be a kernel basis matrix of $\mathbf{A}$, having the form (8), and notice that

[TABLE]

Therefore, any $\boldsymbol{p}\in\mathcal{M}_{0}$, satisfies $\mathbf{D}\log\boldsymbol{p}=\boldsymbol{0}$, and thus, belongs to $RM(\mathbf{A})$. On the other hand, for any $\boldsymbol{p}\in RM({\mathbf{A}})$, both $\bar{\mathbf{D}}\log\boldsymbol{p}=\boldsymbol{0}$ and $\theta_{0}(\boldsymbol{p})=0$ must hold, which immediately implies that $\boldsymbol{p}\in\mathcal{M}_{0}$.

The first claim is proved next.

The relational model $RM({\mathbf{A}})$ is a curved exponential family parameterized by

[TABLE]

If the overall effect is added to $RM({\mathbf{A}})$, the parameter space gets an additional parameter:

[TABLE]

Because $\Theta_{1}$ is the smallest open set in $\mathbb{R}^{J}$ that contains $\Theta$, it parameterizes the minimal regular exponential family containing $RM(\mathbf{A})$. This family is, in fact, $RM(\bar{\mathbf{A}})$. ∎

Example 1**.**

The relational models generated by the matrices

[TABLE]

consist of positive probability distribution which can be written in the following forms:

[TABLE]

where $\beta_{0}$ is the overall effect. The dual representations can be written in the log-linear form, using \mbox{\boldmathd}_{1}=(-1,0,1,-1)^{\prime},\mbox{\boldmathd}_{2}=(-1,1,0,0)^{\prime}\in Ker(\mathbf{A}):

[TABLE]

By Theorem 1, after the overall effect is added, the model specification does not include the non-homogeneous constraint anymore. In terms of the generalized odds ratios:

[TABLE]

The second model may be defined using restrictions only on homogeneous odds ratios, and there is no need to place an explicit restriction on the non-homogeneous odds ratio.

∎

The changes in the model structure after the overall effect is removed are examined next.

A relational model with the overall effect can be reparameterized so that its model matrix has a row of $1$’s, and because of full row rank, this vector is not spanned by the other rows. The implications of the removal of the overall effect will be investigated using a model matrix of this structure, say $\bar{\mathbf{A}}_{1}$. By the removal of the row $\boldsymbol{1}^{\prime}$, one may obtain a different model matrix on the same sample space, but it may happen that there exists a cell $i_{0}$, whose only parameter is the overall effect, and after its removal, the $i_{0}$-th column contains zeros only. In such cases, to have a proper model matrix, such columns, that is such cells, need to be removed. Write $\mathcal{I}_{0}$ for the set of all such cells $i_{0}$, and let $I_{0}=|\mathcal{I}_{0}|$. Then, the reduced model matrix, $\mathbf{A}_{1}$, is obtained from $\bar{\mathbf{A}}_{1}$ after removing the row of $1$’s and deleting the columns which, after this, contain only zeros. This is a model matrix on $\mathcal{I}\setminus\mathcal{I}_{0}$. Without loss of generality, the matrix $\bar{\mathbf{A}}_{1}$ can be written as:

[TABLE]

If the sample spaces of $RM(\bar{\mathbf{A}}_{1})$ and $RM({\mathbf{A}}_{1})$ are the same that is, when $\mathcal{I}_{0}$ is empty, the reduced model is the subset of the original one, consisting of the distributions whose overall effect is zero, see Theorem 2. If the sample space is reduced, the relationship between the kernel basis matrices is described in the next result.

Theorem 3**.**

The following holds:

1. (i)

$dim(Ker(\mathbf{A}_{1}))=dim(Ker(\bar{\mathbf{A}}_{1}))-I_{0}+1$. 2. (ii)

The kernel basis matrix $\mathbf{D}_{1}$ of $\mathbf{A}_{1}$ may be obtained from the kernel basis matrix $\bar{\mathbf{D}}_{1}$ of $\bar{\mathbf{A}}_{1}$ by deleting the the columns in $\mathcal{I}_{0}$ and then leaving out the redundant rows.

Proof.

1. (i)

Because $\bar{\mathbf{A}}_{1}$ is a $J\times I$ matrix of full row rank, $dim(Ker(\bar{\mathbf{A}}_{1}))=I-J$. The linear independence of its rows implies that the rows of $\mathbf{A}_{1}$ are also linearly independent. Therefore, because $\mathbf{A}_{1}$ is a $(J-1)\times(I-I_{0})$ matrix, $dim(Ker(\mathbf{A}_{1}))=I-I_{0}-(J-1)$, which implies the result. 2. (ii)

Let $\boldsymbol{d}_{1},\boldsymbol{d}_{2},\dots,{\boldsymbol{d}}_{I-J}$ be a kernel basis of $\bar{\mathbf{A}}_{1}$. Write

[TABLE]

Then,

[TABLE]

which implies that

[TABLE]

Suppose $\mathbf{A}_{1}$ does not have the overall effect. Notice that each $\boldsymbol{v}_{i}$ has length $I_{0}$, and therefore, one can apply a non-singular linear transformation to the basis vectors $\boldsymbol{d}_{1},\boldsymbol{d}_{2},\dots,{\boldsymbol{d}}_{I-J}$ to reduce them to the form:

[TABLE]

The equations (11) imply that

[TABLE]

[TABLE]

The linear independence of $\boldsymbol{d}_{I_{0}+1},\dots,\boldsymbol{d}_{I-J}$ in $\mathbb{R}^{I}$ entails the linear independence of $\boldsymbol{u}_{I_{0}+1},\dots,\boldsymbol{u}_{I-J}$ in $\mathbb{R}^{I-I_{0}}$. Notice that $\boldsymbol{u}_{1},\dots,\boldsymbol{u}_{I_{0}}$ are jointly linearly independent from $\boldsymbol{u}_{I_{0}+1},\dots,\boldsymbol{u}_{I-J}$, but not necessarily linearly independent from each other. A kernel basis of $\mathbf{A}_{1}$ comprises $I-J-I_{0}+1$ linearly independent vectors in $Ker(\mathbf{A}_{1})$, and, for example, $\boldsymbol{u}_{I_{0}},\boldsymbol{u}_{I_{0}+1},\dots,\boldsymbol{u}_{I-J}$ form such a basis. Therefore, ${\mathbf{D}}_{1}$ can be derived from a kernel basis matrix of $\bar{\mathbf{A}}_{1}$ by removing the columns for $\mathcal{I}_{0}$ and leaving out the $I_{0}-1$ redundant rows.

Suppose $\mathbf{A}_{1}$ has the overall effect; without loss of generality, $\boldsymbol{1}^{\prime}$ is a row of $\mathbf{A}_{1}$. In this case, (11) implies that both $\boldsymbol{1}_{(I-I_{0})}^{\prime}\boldsymbol{u}_{i}=0$ and $\boldsymbol{1}_{I_{0}}^{\prime}\boldsymbol{v}_{i}=0$, for $i=I_{0}+1,\dots,I-J$. Because $\boldsymbol{u}_{i}$’s and $\boldsymbol{v}_{i}$’s vary independently from each other, the linear independence of $\boldsymbol{d}_{1},\boldsymbol{d}_{2},\dots,{\boldsymbol{d}}_{I-J}$ will imply that $(\boldsymbol{u}_{i},\boldsymbol{0})$, for $i=I_{0}+1,\dots,I-J$, are also linearly independent in $\mathbb{R}^{I}$. Consequently, any $I-J-I_{0}+1$ vectors among $\boldsymbol{u}_{i}$’s are linearly independent in $\mathbb{R}^{I-I_{0}}$ and can form a kernel basis of $\mathbf{A}_{1}$. Thus, as in the previous case, ${\mathbf{D}}_{1}$ can be derived from a kernel basis matrix of $\bar{\mathbf{A}}_{1}$ by removing the columns for $\mathcal{I}_{0}$ and leaving out the $I_{0}-1$ redundant rows.

∎

The next two examples illustrate the theorem.

Example 2**.**

Let $RM(\bar{\mathbf{A}}_{1})$ be the relational model generated by

[TABLE]

Here, $\mathcal{I}_{0}=\{5,6\}$. In terms of the generalized odds ratios the model can be written as:

[TABLE]

Remove the row $\boldsymbol{1}^{\prime}$ and the last two columns and consider the reduced matrix:

[TABLE]

The model $RM({\mathbf{A}}_{1})$ does not have the overall effect and can be specified by two generalized odds ratios:

[TABLE]

These odds ratios are defined on the smaller probability space, and may be obtained by removing $p_{5}$ and $p_{6}$, and the redundant odds ratio, from the odds ratio specification of the original model. ∎

Example 3**.**

Consider the relational model $RM(\bar{\mathbf{A}}_{1})$ generated by

[TABLE]

In terms of the generalized odds ratios, the model specification is $p_{1}/p_{3}=1$. Notice that $\bar{\mathbf{A}}_{1}$ is row equivalent to

[TABLE]

Because every $\boldsymbol{d}$ in $Ker(\mathbf{A})$ is orthogonal to $(0,0,0,1)$, its last component has to be zero: $d_{4}=0$. Therefore, in any specification of $RM(\bar{\mathbf{A}}_{1})$ in terms of the generalized odds ratios, $p_{4}$ will not be present. Set

[TABLE]

The model $RM({\mathbf{A}}_{1})$ has the overall effect and can be specified by exactly the same generalized odds ratio as the model $RM(\bar{\mathbf{A}}_{1})$: $p_{1}/p_{3}=1$.

As a further illustration, take

[TABLE]

In this case,

[TABLE]

and $\mathbf{A}_{1}$ is the same as above. In this case, $RM(\bar{\mathbf{A}}_{1})$ is specified by $p_{1}/p_{3}=1,p_{4}/p_{5}=1,p_{4}/p_{6}=1$, and $RM(\bar{\mathbf{A}}_{1})$ is described as previously: $p_{1}/p_{3}=1$. ∎

The polynomial variety $\mathcal{X}_{\bar{\mathbf{A}}_{1}}$ defining the model $RM(\bar{\mathbf{A}}_{1})$ is homogeneous. If the removal of the cells comprising $\mathcal{I}_{0}$ leads to a model without the overall effect, the variety $\mathcal{X}_{\bar{\mathbf{A}}_{1}}$ is dehomogenized, yielding the affine variety $\mathcal{X}_{{\mathbf{A}}_{1}}$ (cf. Cox \BOthers., \APACyear2007).

The converse to this procedure, homogenization of an affine variety, is also studied in algebraic geometry, and is performed by introducing a new variable in such a way that all equations defining the variety become homogeneous. The essence of this procedure is that all probabilities are multiplied by this new variable. This leaves the homogeneous odds ratios unchanged, as the new variable cancels out. The value of a non homogeneous odds ratio becomes, instead of $1$, the reciprocal of the new variable. For example, the odds ratio in Example 2 $p_{2}/(p_{1}p_{4})=1$ becomes $vp_{2}/(vp_{1}vp_{4})=1/v$, where $v$ is the new variable. If now $v$ could be seen as the probability of an additional cell, say $p_{0}$, then this would be a homogeneous odds ratio, $p_{0}p_{2}/(p_{1}p_{4})=1$.

Although a straightforward procedure in algebraic geometry, it does not necessarily have a clear interpretation in statistical inference. Introducing a new variable and a new cell for the purpose of homogenization can be made meaningful in some situations, if the sample space may be extended by one cell, and the new variable is the parameter (probability) of this cell. Homogenization requires this new variable to appear in every cell, too, so the parameter may be seen as the overall effect. The new cell has only the overall effect, thus no feature is present in this cell.

The augmentation of the sample space by an additional cell does make sense, if that cell exists in the population but was not observed because of the design of the data collection procedure, as in Example 4. The additional cell has the overall effect only, thus is a “no feature present” cell.

Example 4**.**

In the study of swimming crabs by Kawamura, Matsuoka, Tajiri, Nishida\BCBL \BBA Hayashi (\APACyear1995), three types of baits were used in traps to catch crabs: fish alone, sugarcane alone, fish-sugarcane combination. The sample space consists of three cells, $\mathcal{I}=\{(0,1),(1,0),(1,1)\}$, and the cell $(0,0)$ is absent by design, because there were no traps without any bait. Under the AS independence, the cell parameter associated with both bait types present is the product of the parameters associated with the other two cells. This is a relational model without the overall effect, generated by the matrix (cf. Klimova \BBA Rudas, \APACyear2015):

[TABLE]

In fact, the overall effect cannot be included in this situation, because it would saturate the model. The affine variety associated with this model can be homogenized by including a new variable. The new variable may only be interpreted as the parameter associated with no bait present, and calls for an additional cell in the sample space (to avoid model the saturation of the model) which may only be interpreted as setting up a trap without any bait, which would also be a plausible research design. The resulting model is generated by $\mathbf{A}_{0}$:

[TABLE]

and indeed, is the model of traditional independence on the complete $2\times 2$ contingency table. ∎

For situations like in Example 4, the AS independence is a natural model, but it also applies to cases, when the “no feature present” situation is logically impossible (like market basket analysis, or records of traffic violations, see Klimova \BOthers. (\APACyear2012); Klimova \BBA Rudas (\APACyear2015), and also the biological example in Section 5), and in such cases, the cell augmentation procedure is not meaningful. There are, however situations, when the existence of the “no feature present” cell is logically not impossible, but the actual existence in the population is dubious.

For a more general discussion of the homogenization of AS independence, let $\boldsymbol{d}_{1},\dots,\boldsymbol{d}_{K}$ be a kernel basis of $\mathbf{A}$, satisfying (8) with $\boldsymbol{d}_{1}^{\prime}\boldsymbol{1}=-1$. The polynomial ideal $\mathscr{I}_{\mathbf{A}}$ associated with the matrix $\mathbf{A}$ is generated by one non-homogeneous polynomial $\boldsymbol{p}^{\boldsymbol{d}_{1}^{+}}-\boldsymbol{p}^{\boldsymbol{d}_{1}^{-}}$, and $K-1$ homogeneous polynomials, $\boldsymbol{p}^{\boldsymbol{d}_{k}^{+}}-\boldsymbol{p}^{\boldsymbol{d}_{k}^{-}}$. Notice that, because $(\boldsymbol{d}_{1}^{+})^{\prime}\boldsymbol{1}-(\boldsymbol{d}_{1}^{-})^{\prime}\boldsymbol{1}$ $=(\boldsymbol{d}_{1}^{+}-\boldsymbol{d}_{1}^{-})^{\prime}\boldsymbol{1}$ = $\boldsymbol{d}_{1}^{\prime}\boldsymbol{1}$, the difference in the degrees of the monomials $p^{\boldsymbol{d}_{1}^{+}}$ and $p^{\boldsymbol{d}_{1}^{-}}$ is $-1$. Therefore, the polynomial $\boldsymbol{p}^{\boldsymbol{d}_{1}^{+}}-\boldsymbol{p}^{\boldsymbol{d}_{1}^{-}}$ can be homogenized by multiplying the first monomial by one additional variable, say $p_{0}$:

[TABLE]

The polynomial ideal generated by

[TABLE]

and the corresponding variety are homogeneous, and can be described by the matrix of size $(J+1)\times(I+1)$ of the following structure:

[TABLE]

Here, $\boldsymbol{1}^{\prime}_{I}$ is the row of $1$’s of length $I$, and $\boldsymbol{0}_{J}$ is the column of zeros of length $J$.

In fact, the homogeneous variety $\mathcal{X}_{{\mathbf{A}}_{0}}$ is the projective closure of the affine variety $\mathcal{X}_{\bar{\mathbf{A}}_{0}}$ (Cox \BOthers., \APACyear2007). The latter can be obtained from the former by dehomogenization via setting $p_{0}=1$.

The homogenization of the model of AS independence for three features is discussed next.

Example 5**.**

Consider the model of AS independence for three attributes, $A$, $B$, and $C$, described in Klimova \BOthers. (\APACyear2012).

[TABLE]

Here $p_{ijk}=\mathbb{P}(A=i,B=j,C=k)$ for $i,j,k\in\{0,1\}$, where the combination $(0,0,0)$ does not exist, and $\sum_{ijk}p_{ijk}=1$. The equations (12) specify the relational model generated by

[TABLE]

Consider the following kernel basis matrix which is of the form (8):

[TABLE]

The corresponding polynomial ideal is:

[TABLE]

The generating set of $\mathscr{I}_{\mathbf{A}}$ includes at least one non-homogeneous polynomial, due to $\boldsymbol{d}_{1}$, and can be homogenized by introducing a new variable, say $p_{000}$. The resulting ideal,

[TABLE]

is homogeneous, and its zero set

[TABLE]

where

[TABLE]

is thus a homogeneous variety. The relational model $RM(\mathbf{A}_{0})$ is defined on a larger sample space, namely $\mathcal{I}\cup(0,0,0)$. The model has the overall effect and is the following set of distributions:

[TABLE]

The rows of $\mathbf{A}_{0}$ are the indicators of the cylinder sets of the total (the row of $1$’s), and of the $A$, $B$, and $C$ marginals. Therefore, the relational model $RM(\mathbf{A}_{0})$ is the traditional model of mutual independence. ∎

The next theorem states in general what was seen in the example. Let $X_{1},\dots,X_{T}$ be the random variables taking values in $\{0,1\}$. Write $\mathcal{I}^{0}$ for the Cartesian product of their ranges, and let $\mathcal{I}=\mathcal{I}^{0}\setminus(0,\dots,0)$.

Theorem 4**.**

Let $RM(\mathbf{A})$ be the model of AS independence of $X_{1},\dots,X_{T}$ on the sample space $\mathcal{I}$. The interior of the projective closure of this model is the log-linear model of mutual independence of $X_{1},\dots,X_{T}$ on the sample space $\mathcal{I}^{0}$.

Proof.

Let $\mathbf{A}$ be the model matrix for the AS independence:

[TABLE]

The number of columns of $\mathbf{A}$ is equal to the number of cells in the sample space $\mathcal{I}$, $I=2^{T}-1$. The model $RM(\mathbf{A})$ is the intersection of the polynomial variety $\mathcal{X}_{\mathbf{A}}$ and the interior of the simplex $\Delta_{I-1}$. The variety $\mathcal{X}_{\mathbf{A}}$ is non-homogeneous, because among its generators there is at least one non-homogeneous polynomial. In order to obtain the projective closure of $\mathcal{X}_{\mathbf{A}}$ (cf. Cox \BOthers., \APACyear2007), include the “no feature present” cell, indexed by [math], to the sample space, choose a Gröbner basis of the ideal $\mathscr{I}_{\mathbf{A}}$, and homogenize all non-homogeneous polynomials in this basis using the cell probability $p_{0}$. Because the projective closure of $\mathcal{X}_{\mathbf{A}}$ is the minimal homogeneous variety in the projective space whose dehomogenization is $\mathcal{X}_{\mathbf{A}}$ (Cox \BOthers., \APACyear2007), Theorem 3(ii) implies that this closure can be described using the matrix

[TABLE]

Each distribution in $RM(\mathbf{A})$ has the multiplicative structure prescribed by $\mathbf{A}$ (Klimova \BBA Rudas, \APACyear2016), and during the homogenization, is mapped in a positive distribution in $\mathcal{X}_{\mathbf{A}_{0}}$. Because all strictly positive distributions in $\mathcal{X}_{\mathbf{A}_{0}}$ have the multiplicative structure prescribed by $\mathbf{A}_{0}$, they comprise the relational model $RM(\mathbf{A}_{0})$. This matrix describes the model of mutual independence between $X_{1},\dots,X_{T}$ in the effect coding, and the proof is complete. ∎

The homogenization (in the language of algebraic geometry) or regularization (in the language of the exponential families) leads to a simpler structure, which allows a simpler calculation of the MLE. However, the additional cell was not observed in these cases, and assuming its frequency is zero is ungrounded and may lead to wrong inference.

The framework developed here may also be used to define context specific independence, so that in one context conditional independence holds, in another one, AS independence does. To illustrate, let $X_{1}$, $X_{2}$, $X_{3}$ be random variables taking values in $\{0,1\}$. Assume that the $(0,0,0)$ outcome is impossible, so the sample space can be expressed as:

[TABLE]

Let $\boldsymbol{p}=(p_{001},p_{010},p_{011},p_{100},p_{101},p_{110},p_{111})$, and consider the relational model without the overall effect generated by

[TABLE]

The kernel basis matrix is equal to:

[TABLE]

and thus, the model can be specified in terms of the following two generalized odds ratios:

[TABLE]

The second constraint expresses the (conventional) context-specific independence of $X_{1}$ and $X_{2}$ given $X_{3}=1$. The first odds ratio is non-homogeneous, and the corresponding constraint may be seen as the context-specific AS-independence of of $X_{1}$ and $X_{2}$ given $X_{3}=0$.

4 ML estimation with and without the overall effect

The properties of the ML estimates under relational models, discussed in detail in Klimova \BOthers. (\APACyear2012) and Klimova \BBA Rudas (\APACyear2016), are summarized here in the language of the linear and multiplicative families defined by the model matrix and its kernel basis matrix. The conditions of existence of the MLE are reviewed first.

Let $\boldsymbol{a}_{1},\dots,\boldsymbol{a}_{|\mathcal{I}|}$ denote the columns of $\mathbf{A}$, and let $C_{\mathbf{A}}=\{\boldsymbol{t}\in\mathbb{R}^{J}_{\geq 0}:\,\,\exists\boldsymbol{p}\in\mathbb{R}^{|\mathcal{I}|}_{\geq 0}\quad\boldsymbol{t}=\mathbf{A}\boldsymbol{p}\}$ be the polyhedral cone whose relative interior comprises such $\boldsymbol{t}\in\mathbb{R}^{J}_{>0}$, for which there exists a $\boldsymbol{p}>\boldsymbol{0}$ that satisfies $\boldsymbol{t}=\mathbf{A}\boldsymbol{p}$. A set of indices $F=\{i_{1},i_{2},\dots,i_{f}\}$ is called facial if the columns $\boldsymbol{a}_{i_{1}},\boldsymbol{a}_{i_{2}},\dots,\boldsymbol{a}_{i_{f}}$ are affinely independent and span a proper face of $C_{\mathbf{A}}$ (cf. Grünbaum, \APACyear2003; Geiger, Meek\BCBL \BBA Sturmfels, \APACyear2006; Fienberg \BBA Rinaldo, \APACyear2012). It can be shown that a set $F$ is facial if and only if there exists a $\boldsymbol{c}\in\mathbb{R}^{J}$, such that $\boldsymbol{c}^{\prime}\boldsymbol{a}_{i}=0$ for every $i\in F$ and $\boldsymbol{c}^{\prime}\boldsymbol{a}_{i}>0$ for every $i\notin F$.

Let $\boldsymbol{q}\in\mathcal{P}$ and let $\mathcal{K}$ be the set of $\kappa>0$, such that, for a fixed $\kappa$, the linear family

[TABLE]

is not empty, and let \mathcal{F}(\mathbf{A},\mbox{\boldmathq})=\bigcup_{\mathcal{K}}\mathcal{F}(\mathbf{A},\mbox{\boldmathq},\kappa). For each $\kappa>0$, the linear family \mathcal{F}(\mathbf{A},\mbox{\boldmathq},\kappa) is a polyhedron in the cone $C_{\mathbf{A}}$.

Theorem 5**.**

(Klimova \BBA Rudas, \APACyear2016)* Let $RM(\mathbf{A})$ be a relational model, with or without the overall effect, and let $\boldsymbol{q}$ be the observed distribution.*

The MLE $\hat{p}_{\boldsymbol{q}}$ given $\boldsymbol{q}$ exists if only:

1. (i)

$supp(\boldsymbol{q})=\mathcal{I}$, or 2. (ii)

$supp(\boldsymbol{q})\subsetneq\mathcal{I}$* and, for all facial sets $F$ of $\mathbf{A}$, $supp(\boldsymbol{q})\not\subseteq F$. *

In either case, \hat{\mbox{\boldmathp}}_{q}=\mathcal{F}(\mathbf{A},\mbox{\boldmathq})\cap int(\mathcal{X}_{\mathbf{A}}), and there exists a unique constant $\gamma_{q}>0$, also depending on $\mathbf{A}$, such that:

[TABLE] 2. 2.

The MLE under the extended model \,\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-2.5ptRM\kern-1.00006pt}}}(\mathbf{A}), defined in (7), always exists and is the unique point of $\mathcal{X}_{\mathbf{A}}$ which satisfies:

[TABLE]

The statements follow from Theorem 4.1 in Klimova \BBA Rudas (\APACyear2016) and Corollary 4.2 in Klimova \BOthers. (\APACyear2012), and the proof is thus omitted. The constant $\gamma_{q}$, called the adjustment factor, is the ratio between the subset sums of the MLE, \mathbf{A}\hat{\mbox{\boldmathp}}_{q}, and the subset sums of the observed distribution, \mathbf{A}\mbox{\boldmathq}. If the overall effect is present in the model, $\gamma_{q}=1$ for all $\boldsymbol{q}$.

Let $\mathbf{A}$ be a model matrix whose row space does not contain $\boldsymbol{1}^{\prime}$, and let $\bar{\mathbf{A}}$ be the matrix obtained by augmenting $\mathbf{A}$ with the row $\boldsymbol{1}^{\prime}$. It will be shown in the proof of the next theorem that every facial set of $\mathbf{A}$ is facial for $\bar{\mathbf{A}}$. If the observed $\boldsymbol{q}$ is positive, the MLEs \hat{\mbox{\boldmathp}}_{q} and \bar{\mbox{\boldmathp}}_{q} under the models $RM({\mathbf{A}})$ and $RM(\bar{\mathbf{A}})$, respectively, exist. However, as implied by the relationship between the facial sets of ${\mathbf{A}}$ and $\bar{\mathbf{A}}$, if $\boldsymbol{q}$ has some zeros, the MLE may exist under $RM({\mathbf{A}})$, but not under $RM(\bar{\mathbf{A}})$, or neither of the MLEs exist.

Theorem 6**.**

Let $\mathbf{A}$ be a model matrix whose row space does not contain $\boldsymbol{1}^{\prime}$, and let $\bar{\mathbf{A}}$ be the matrix obtained by augmenting $\mathbf{A}$ with the row $\boldsymbol{1}^{\prime}$. Let $\boldsymbol{q}$ be the observed distribution. If, given $\boldsymbol{q}$, the MLE under $RM(\bar{\mathbf{A}})$ exists, so does the MLE under $RM({\mathbf{A}})$.

Proof.

If $\boldsymbol{q}>\boldsymbol{0}$, both MLEs exists.

Assume that $\boldsymbol{q}$ has some zeros, that is, $supp(\boldsymbol{q})\subsetneq\mathcal{I}$, and that the MLE under $RM(\bar{\mathbf{A}})$ exists. It will be shown next that for any facial set $F$ of $\mathbf{A}$, $supp(\boldsymbol{q})\not\subseteq F$.

The proof is by contradiction. Let $F_{0}$ be a facial set of $\mathbf{A}$, such that $supp(\boldsymbol{q})\subset F_{0}$. Therefore, there exists a $\boldsymbol{c}\in\mathbb{R}^{J}$, such that $\boldsymbol{c}^{\prime}\boldsymbol{a}_{i}=0$ for every $i\in F$ and $\boldsymbol{c}^{\prime}\boldsymbol{a}_{i}>0$ for every $i\notin F$.

Denote by $\bar{\boldsymbol{a}}_{1},\dots,\bar{\boldsymbol{a}}_{I}$ the columns of $\bar{\mathbf{A}}$. By construction, $\bar{\boldsymbol{a}}_{i}=(1,\boldsymbol{a}_{i}^{\prime})^{\prime}$, $\,i=1,\dots,I$. Let $\bar{\boldsymbol{c}}=(0,\boldsymbol{c}^{\prime})^{\prime}$. Then,

[TABLE]

and thus, $F_{0}$ is a facial set of $\bar{\mathbf{A}}$. Because $supp(\boldsymbol{q})\subset F_{0}$, the MLE under $RM(\bar{\mathbf{A}})$, given $\boldsymbol{q}$, does not exist, which contradicts the initial assumption. This completes the proof. ∎

Example 5 (revisited): Let $\boldsymbol{q}_{1}=(0,0,0,0,0,0,1)^{\prime}$ be the observed distribution. Because $supp(\boldsymbol{q}_{1})=\{7\}$ is not a subset of any facial sets of $\mathbf{A}$, the MLE exists:

[TABLE]

with $\hat{\gamma}_{q}=2-\sqrt[3]{4}$.

On the other hand, the set of indices $F=\{1,4,5,7\}$ is facial for $\bar{\mathbf{A}}$, and $supp(\boldsymbol{q}_{1})\subsetneq F$. In this case, the MLE exists only in the extended model \,\hbox{\vbox{\hrule height=0.5pt\kern 2.15277pt\hbox{\kern-2.5ptRM\kern-1.00006pt}}}(\bar{\mathbf{A}}), and is equal to $\boldsymbol{q}_{1}$ itself.

Let $\boldsymbol{q}_{2}=(1,0,0,0,0,0,0)^{\prime}$. Because $supp(\boldsymbol{q}_{2})=\{1\}$ is a subset of a facial set of $\mathbf{A}$ and of a facial set of $\bar{\mathbf{A}}$, the MLEs exist only in the corresponding extended models. ∎

Further properties of the adjustment factor, including its geometrical meaning, are described next, relying on the following result:

Theorem 7**.**

Let $\mathbf{A}$ be a model matrix whose row space does not contain $\boldsymbol{1}^{\prime}$, and let $\bar{\mathbf{A}}$ be the matrix obtained by augmenting $\mathbf{A}$ with the row $\boldsymbol{1}^{\prime}$. For any \mbox{\boldmathr}_{1},\>\mbox{\boldmathr}_{2}\in\mathcal{P}, \mbox{\boldmathr}_{1}\neq\mbox{\boldmathr}_{2}, the following holds:

2. (i)

The MLEs under $RM(\mathbf{A})$, given they exist, are equal if and only if the subset sums entailed by $\mathbf{A}$ are proportional:

[TABLE]

and the adjustment factors in the MLE satisfy: $\kappa\hat{\gamma}_{r_{1}}$ = $\hat{\gamma}_{r_{2}}$. 2. 2.

The MLEs under $RM(\bar{\mathbf{A}})$, given they exist, are equal if and only if the subset sums entailed by $\mathbf{A}$ coincide:

[TABLE]

The statements are a reformulation of Theorem 4.4 in Klimova \BOthers. (\APACyear2012), and no proofs are provided here. The relationship between the adjustment factors is obvious.

The theorem implies that \mathcal{F}(\mathbf{A},\mbox{\boldmathq}) is an equivalence class in $\mathcal{P}$, in the sense that, for any \mbox{\boldmathr}\in\mathcal{F}(\mathbf{A},\mbox{\boldmathq}), the MLE under $RM(\mathbf{A})$ satisfies \hat{\mbox{\boldmathp}}_{r} = \hat{\mbox{\boldmathp}}_{q}. Each sub-family \mathcal{F}(\mathbf{A},\mbox{\boldmathq},\kappa) is characterized by its unique adjustment factor under $RM(\mathbf{A})$. That is, for every $\boldsymbol{r}_{1},\boldsymbol{r}_{2}\in\mathcal{F}(\mathbf{A},\boldsymbol{q},\kappa)$, $\,\,\boldsymbol{r}_{1}\neq\boldsymbol{r}_{2}$,

[TABLE]

In addition, $\,\,\bar{\boldsymbol{p}}_{r_{1}}=\bar{\boldsymbol{p}}_{r_{2}}$ for any $\boldsymbol{r}_{1},\boldsymbol{r}_{2}\in\mathcal{F}(\mathbf{A},\boldsymbol{q},\kappa)$, and therefore, for a fixed $\kappa$, \mathcal{F}(\mathbf{A},\mbox{\boldmathq},\kappa) is an equivalence class under $RM(\bar{\mathbf{A}})$.

From a geometrical point of view, \mathcal{F}(\mathbf{A},\mbox{\boldmathq}) is a polyhedron which decomposes into polyhedra \mathcal{F}(\mathbf{A},\mbox{\boldmathq},\kappa), with $\kappa>0$; clearly, \mbox{\boldmathq}\in\mathcal{F}(\mathbf{A},\mbox{\boldmathq},1). The MLE under $RM(\bar{\mathbf{A}})$ given \mbox{\boldmathr}\in\mathcal{F}(\mathbf{A},\mbox{\boldmathq},\kappa) is the unique point common to the polyhedron \mathcal{F}(\mathbf{A},\mbox{\boldmathq},\kappa) and the variety $\mathcal{X}_{\bar{\mathbf{A}}}$. Among the feasible values of $\kappa$ there exists a unique one, say $\hat{\kappa}$, such that the MLE \bar{\mbox{\boldmathp}}_{r}, \forall\mbox{\boldmathr}\in\mathcal{F}(\mathbf{A},\mbox{\boldmathq},\hat{\kappa}), coincides with the MLE of $q$ under $RM(\mathbf{A})$, \hat{\mbox{\boldmathp}}_{q}. This happens when $\hat{\gamma}_{r}=1$ so that, from (ii) in Theorem 7, $\hat{\kappa}$ = $\hat{\gamma}_{q}$. This latter point, \hat{\mbox{\boldmathp}}_{q}, is the intersection between \mathcal{F}(\mathbf{A},\mbox{\boldmathq}) and the non-homogeneous variety $\mathcal{X}_{{\mathbf{A}}}$. This specific value of the adjustment factor $\gamma_{q}=\hat{\kappa}$, is the adjustment factor of the MLE under $RM(\mathbf{A})$ given $q$. An illustration is given next.

Relational models for probabilities without the overall effect are curved exponential families, and the computation of the MLE under such models is not straightforward. An extension of the iterative proportional fitting procedure, G-IPF, that can be used for both models with and models without the overall effect was proposed in Klimova \BBA Rudas (\APACyear2015) and is implemented in Klimova \BBA Rudas (\APACyear2014). Alternatively, the MLEs can be computed, for instance, using the Newton-Raphson algorithm or the algorithm of Evans \BBA Forcina (\APACyear2013). One of the algorithms described in Forcina (\APACyear2017) gave an idea of a possible modification of G-IPF. A brief description of the original and modified versions of G-IPF is given below:

[TABLE]

Theorem 8**.**

If $\boldsymbol{q}>\boldsymbol{0}$, the G-IPFm algorithm converges, and its limit is equal to $\hat{\boldsymbol{p}}_{q}$, the ML estimate of $\boldsymbol{p}$ under $RM(\mathbf{A})$.

Proof.

The convergence of one iteration of G-IPFm, when $\gamma$ is fixed, can be proved similarly to Theorem 3.2 in Klimova \BBA Rudas (\APACyear2015). The limit is positive, $\tilde{\boldsymbol{p}}_{\gamma}>\boldsymbol{0}$, and thus, by Lemma 1 in Forcina (\APACyear2017)f(\gamma)=\boldsymbol{d}_{1}^{\prime}\log{\tilde{\mbox{\boldmathp}}}_{\gamma} is a strictly increasing and differentiable function of $\gamma$. So, one can update $\gamma$, until for some ${\gamma}_{q}$ the G-IPFm limit satisfies: $f(\gamma_{q})=\boldsymbol{d}_{1}\log\tilde{\boldsymbol{p}}_{{\gamma}_{q}}=0$. Because, in this case,

[TABLE]

the uniqueness of the MLE implies that \tilde{\mbox{\boldmathp}}_{\gamma_{q}}=\hat{\mbox{\boldmathp}}_{q} and ${\gamma}_{q}=\hat{\gamma}_{q}$. ∎

The original G-IPF can be used whether or not $\boldsymbol{q}$ has some zeros, and it computes a sequence whose elements are the unique intersections of the variety $\mathcal{X}_{\mathbf{A}}$ and each of the polyhedra defined by \mathbf{A}\tilde{\mbox{\boldmath\tau}}=\gamma\mathbf{A}\mbox{\boldmathq} for different $\gamma$. This sequence converges, and its limit belongs to the hyperplane $\boldsymbol{1}^{\prime}\boldsymbol{\tau}=1$ (Klimova \BBA Rudas, \APACyear2016). G-IPFm produces a sequence whose elements are the unique intersections of the interior of the homogeneous variety $\mathcal{X}_{\bar{\mathbf{A}}}$ and each of the polyhedra \mathcal{F}(\mathbf{A},\mbox{\boldmathq},\gamma). The limit of this sequence belongs to the interior of the non-homogeneous variety $\mathcal{X}_{\mathbf{A}}$. To ensure the existence, differentiability, and monotonicity of $f(\gamma)$, described above, the G-IPFm algorithm should be applied only when $\boldsymbol{q}>\boldsymbol{0}$. If $\boldsymbol{q}$ has some zero components, the positive MLE \hat{\mbox{\boldmathp}}_{q} may still exist, see Theorem 5(ii). However, for some $\boldsymbol{q}$, because, in general, the matrices $\mathbf{A}$ and $\bar{\mathbf{A}}$ have different facial sets, no strictly positive $\boldsymbol{p}_{\gamma}$ would satisfy \bar{\mathbf{A}}\boldsymbol{p}_{\gamma}=\left(\begin{array}[]{c}1\\ \gamma\mathbf{A}\boldsymbol{q}\end{array}\right).

Some limitations and advantages of using the generalized IPF were addressed in Klimova \BBA Rudas (\APACyear2015), Section 2. In particular, while the assumption of the model matrix to be of full row rank can be relaxed for G-IPF, it is one of the major assumptions for the Newton-Raphson and the Fisher scoring algorithms. The algorithms proposed in Forcina (\APACyear2017) also require the model matrix to be of full row rank, and their convergence relies on the positivity of the observed distribution.

5 Loss of potentials in hematopoiesis

Hematopoietic stem cells (HSC) are able to become progenitors that, in turn, may develop into mature blood cells. Understanding the process of forming mature blood cells, called hematopoiesis, is one of the most important aims of cell biology, as it may help to develop new cancer treatments. The HSC progenitors can proliferate (produce cells of the same type) or differentiate (produce cells of different types). Multiple experiments suggested that HSC progenitors are multipotent cells and differentiate by losing one of the potentials. While the mature blood cells are unipotent, they do not proliferate or differentiate The differentiation is believed to be a hierarchical process, with HSC progenitors and mature blood cells at the highest and the lowest levels, respectively.

The models discussed below apply to the steady-state of hematopoiesis, under the assumption that cells neither proliferate nor die and can undergo only first phase of differentiation. Various hierarchical models for differentiation have been proposed (cf. Kawamoto, Wada\BCBL \BBA Katsura, \APACyear2010; Ye, Huang\BCBL \BBA Guo, \APACyear2017). The equal loss of potentials (ELP) model was introduced in Perié \BOthers. (\APACyear2014), and is described next. Denote by $MDB$ the three-potential HSC progenitor of the $M$, $D$, and $B$ mature blood cell types. During the first phase of differentiation, an $MDB$ progenitor can differentiate by losing either one or two potentials at the same time, and thus produce a cell of one of the six types: $M$, $D$, $B$, $MD$, $MB$, $DB$.

Let $p$ be the vector of probabilities of losing the corresponding potentials from $MDB$:

[TABLE]

For example, $p_{*DB}$ is the probability of losing the $M$ potential from $MDB$, $p_{M*B}$ is the probability of loosing the $D$ potential from $MDB$, and $p_{**B}$ is the probability of losing the $M$ and $D$ potentials from $MBD$ at the same time, and so on. The ELP model assumes that “the probability to lose two potentials at the same time is the product of the probability of losing each of the potentials” (see Caption to Fig 3A, Perié \BOthers., \APACyear2014):

[TABLE]

The model specified by (19) is the relational model generated by the matrix

[TABLE]

or, in a parametric form,

[TABLE]

where, using the notation in Perié \BOthers. (\APACyear2014), $\alpha_{M},\alpha_{D},\alpha_{B}$ are the parameters associated with the loss of the corresponding potential from $MDB$. It can be easily verified that the relational model generated by (20) does not have the overall effect, so the normalization has to be added as a separate condition:

[TABLE]

Perié \BOthers. (\APACyear2014) define the ELP model in the following parametric form:

[TABLE]

That is, the authors rescaled the loss probabilities to force them sum to $1$. In fact, (5) is also a relational model; it is generated by

[TABLE]

and can be obtained by adding the overall effect to the model defined by (20). Because, the original model does not have the overall effect, adding a row of $1$’s changed this model. One can check by substitution that the probabilities in (5) do not satisfy the multiplicative constraints (19). The estimates of the probabilities of loss of potentials from the $MDB$ cells shown in Figure 3B of Perié \BOthers. (\APACyear2014). In the notation used here,

[TABLE]

These probabilities sum to $1$, but also do not satisfy (19).

Acknowledgments

The authors wish to thank Antonio Forcina for his thought-provoking discussions, Ingmar Glauche and Christoph Baldow for their help with understanding the main concepts of hematopoiesis, and Wicher Bergsma. The second author is also a Recurrent Visiting Professor at the Central European University and the moral support received is acknowledged.