A biclustering algorithm for binary matrices based on penalized Bernoulli likelihood

Abstract

We propose a new biclustering method for binary data matrices using the maximum penalized Bernoulli likelihood estimation. Our method applies a multi-layer model defined on the logits of the success probabilities, where each layer represents a simple bicluster structure and the combination of multiple layers is able to reveal complicated, multiple biclusters. The method allows for non-pure biclusters, and can simultaneously identify the 1-prevalent blocks and 0-prevalent blocks. A computationally efficient algorithm is developed and guidelines are provided for specifying the tuning parameters, including initial values of model parameters, the number of layers, and the penalty parameters. Missing-data imputation can be handled in the EM framework. The method is tested using synthetic and real datasets and shows good performance.

Appendix

1.1 Derivations of (4) and (5)

Lemma 1

The function π(x){1−π(x)} is decreasing in x≥0 where π(x)={1+exp(−x)}⁻¹.

Proof

The first derivative of the underlying function is π′(x){1−π(x)}−π(x)π′(x)=π′(x){1−2π(x)}=π(x){1−π(x)}{1−2π(x)}. Since 1/2≤π(x)≤1 for x≥0, the derivative is negative and therefore the function is decreasing. □

Lemma 2

The function $r(x)=\log\pi(\sqrt{x})-\sqrt{x}/2$ is convex.

Proof

The second derivative of r(x) is

$$r''(x) = \frac{1}{4x} \biggl[ \frac{2\pi(\sqrt{x})-1}{2\sqrt{x}} - \pi(\sqrt{x})\bigl\{1-\pi(\sqrt{x})\bigr\} \biggr]. $$

Note that

with $\xi\in(-\sqrt{x},\sqrt{x})$ from the mean value theorem. From $\xi<\sqrt{x}$ and Lemma 1, the second derivative of r(x) is positive. This completes the proof of Lemma 2. □

From the convexity of r(x), we get r(x)≥r(y)+r′(y)(x−y) at any y, so that

or

By replacing $\sqrt{x}$ and $\sqrt{y}$ with x and y respectively, we obtain

The coefficient of the quadratic term is bounded above by 1/8 because

for ξ∈(−y,y) by the mean value theorem. At y=0, the coefficient of the quadratic term is not defined properly. In this case, we use the limit when y approaches zero. By L’hopital’s theorem we get

$$\lim_{y\rightarrow0} \frac{2\pi(y)-1}{4y} = \lim_{y\rightarrow 0}\frac {2\pi'(y)}{4} = \lim_{y\rightarrow0}\frac{\pi(y)\{1-\pi(y)\}}{2} = \frac{1}{8}. $$

Now, by completing squares around x, we get the upper bound as

$$-\log\pi(y) + 2\bigl\{1-\pi(y)\bigr\}^2 + \frac{1}{8}\bigl[x- \bigl\{y+4\bigl(1-\pi(y)\bigr)\bigr\}\bigr]^2. $$

Replacing x and y by q _ij θ _ij and $q_{ij}\theta_{ij}^{o}$ respectively, we obtain the upper bound in (4). In fact, the first two terms in (4) are obtained straightforwardly after x and y are replaced by q _ij θ _ij and $q_{ij}\theta _{ij}^{o}$ respectively. After the replacement, the third term of the above displayed formula becomes

In the above, we used $q_{ij}^{2}=1$ since q _ij takes values −1 or 1 only. Now, if we define $x_{ij}=\theta_{ij}^{o}+4q_{ij} \{1-\pi (q_{ij}\theta_{ij}^{o}) \}$, the desired result is obtained.

Applying (4) to (1), we obtain the upper bound of the penalized Bernoulli log likelihood

where $g(\mbox{\boldmath$\varXi$}|\mbox{\boldmath$\varXi$}^{o})$ is in the equation (5) of the manuscript, and C is the constant term $-\sum_{i=1}^{n}\sum _{j=1}^{p}\log (q_{ij}\theta_{ij}^{o}) + 2\sum_{i=1}^{n}\sum_{j=1}^{p}\{1-\pi (q_{ij}\theta _{ij}^{o})\}^{2}$, which does not involve θ _ij . Thus, g can be used as the surrogate function for the upper-bound optimization.