1]Department of Astronomy and Bosscha Observatory, FMIPA, Institut Teknologi
Bandung, Jalan Ganesha 10, Bandung 40132, Indonesia 2] SYRTE, Observatoire de
Paris, PSL Research University, CNRS UMR 8630, Paris, France ] Email:
<EMAIL_ADDRESS>
# Study of the relation between luni-solar periodicities and earthquake events

Ibnu NURUL HUDA  Jean SOUCHAY  [  [  [

(2 December 2020)

###### Abstract

We revisit the work of Kilston and Knopoff (1983) to study the correlation of
earthquakes with the main lunisolar tidal components in the limited zone of
the Southern California region, with a considerably bigger amount of data. By
adopting the same methodology as these authors, we confirm their results, i.e.
that strong earthquakes with M $\geq 6$ show a statistically significant
correlation with half day and precession cycles, with respective periods 12
hr. and 18.6 y., while weaker earthquakes show a correlation with the annual
cycle. Moreover we extended the study to the analysis of bi-modal phase to
phase diagrams, leading to a very high probability of a lack of earthquake
occurrences when both phases of the 12 hr. and 18.6 y. cycles are inside the
[0 - 0.5] interval. Furthermore, we demonstrate that the statistical method is
not really appropriate to study the correlation with the semi-lunar day and
the fortnightly tidal components.

###### keywords:

Earthquake interaction – Tides and planetary waves – Statistical test

##  1 Introduction

The relation between lunisolar cycle and earthquake occurrences has been a
long-standing debate for more than one century. Indeed the Earth tide produced
by the Moon and the Sun is subject to very well modelized oscillations and the
corresponding tidal stress is superimposed on tectonic stress. Although the
amplitude of the tidal stress changes is generally considerably smaller than
the stress drop related to earthquakes, typically of the order of $10^{3}$
Pa., its frequency is generally much larger. Therefore when the stress level
is critical with respect to the release of an earthquake, the rapid additional
stress change due to the Earth tide may trigger an earthquake.

Based on this principle, there have been numerous searches for the relation
between Earth tides and earthquake occurrences. Some authors concluded that
there is a correlation between these two phenomena (see, e.g., Ryall,
VanWormer, and Jones, 1968; Tamrazyan, 1974; Young et al., 1978; Nikolaev,
1995; Saltykov, Ivanov, and Kugaenko, 2004) whereas others reject this
conclusion (see, e.g., Knopoff, 1964; Shlien, 1972; Shudde and Barr, 1977). In
general, the correlation of earthquakes with lunisolar periodicities is
studied in two ways. First, several authors have devoted their efforts to the
correlation on a global scale (see, e.g., Tanaka, Ohtake, and Sato, 2002;
Métivier et al., 2009). Meanwhile, since the nature of the correlations seems
to be strongly dependent on the faulting mechanism involved, other authors
conducted the study by considering a limited zone rather than a global one.
Such local investigations were done for instance for the Southern California
area (Kilston and Knopoff, 1983), the Italy territory (Palumbo, 1986), the
Kamchatka area (Shirokov, 1983), the Vrancea area in Romania (Souchay and
Stavinschi, 1997; Stavinschi and Souchay, 2003; Cadicheanu, Van Ruymbeke, and
Zhu, 2007), Taiwan (Chen et al., 2012), Hellenic Arc (Vergos, Arabelos, and
Contadakis, 2015), and Colombia (Moncayo, Zuluaga, and Monsalve, 2019).

Here we revisit the work of Kilston and Knopoff (1983), hereafter named as K &
K, to study the earthquake correlation with lunisolar cycles in the limited
zone of the Southern California. Following this last paper, the reason to
search such a correlation stays in geotectonic considerations: as it is the
case for the San Andreas fault, the majority of faults in this region are
aligned along an NW-SE direction, and their characteristics are right-lateral
strike-slip. As a result, extensional tidal stresses in the east-west
direction increase the shear stress at the origin of the faulting. Note that
Earth tides are the largest known sources of oscillatory stresses with
relatively short periods, as it is asserted by K & K. Since east-west
extensional stress is maximized for the Moon and the Sun respectively at the
Moonrise/set and the Sunrise/set we can consider the intuitive hypothesis that
earthquakes occur preferentially during these instants, as those corresponding
to maxima of other leading tidal components. It should be all the more the
case when both bodies, Moon and Sun are close to the configuration of
syzygies, that is to say near the new Moon and the full Moon, when their
gravitational action is combined. By using the same statistical approach as K
& K, we take advantage of the better quality of the seismologic data available
as well as new data gathered during roughly four more decades after K & K
publication.

This paper is constructed as follows. In Section 2, we explain in detail both
the contents of the seismologic data used and the statistical method based on
$\chi^{2}$ analysis. The phase distribution of the earthquake events relative
to the main lunisolar tidal cycles is analyzed in Section 3. Section 4
presents the estimation of correlations between Earth tides components and
earthquake occurrences through the same statistical test as used by K & K.
Section 5 discusses the sensitivity of the results with respect to the period
of the cycle considered. Finally, the conclusion is given in Section 6.

##  2 Data and method

We use the earthquake data from United States Geological Survey
(https://earthquake.usgs.gov/), by using, in particular, the catalogues CI
(Center, 2013) and USHIS (Stover and Coffman, 1993), with the data starting
from 1800 to 2020. We limit our investigations to the region inside the
latitude interval from $33.0^{\circ}$ to $36.0^{\circ}$ and the longitude
interval westward from $-115.5^{\circ}$ to $-120.3^{\circ}$. To minimize the
possible contamination of our data set by aftershocks, we limit the study by
considering only the earthquakes with magnitude M $\geq 5.3$ and eliminate
earthquake aftershocks that occur close to strong earthquakes M $\geq 6$. Thus
by using this restriction, we get 68 isolated earthquake events. 32 of them
have a magnitude M $\geq 6$ and the remaining 36 ones range in the interval
$5.3\leq M<6$. These numbers are almost twice bigger than in the data used by
K & K. Tables 1 and 2 display respectively the earthquake events found in both
magnitude intervals.

Table 1: List of earthquake events with M $\geq 6$. The symbol * signifies that the event was taken into account by K & K. Date | hr:min |  lat (∘) |  long (∘) | mag
---|---|---|---|---
| (UTC) |  |  |
2019-07-06 | 03:19 | 35.8 | -117.6 | 7.1
2010-07-07 | 23:53 | 33.4 | -116.5 | 6.0
1999-10-16 | 09:46 | 34.6 | -116.3 | 6.6
1994-01-17 | 12:30 | 34.2 | -118.5 | 6.4
1992-06-28 | 11:57 | 34.2 | -116.4 | 7.3
1987-11-24 | 13:15 | 33.0 | -115.8 | 6.0
1971-02-09 * | 14:00 | 34.4 | -118.4 | 6.4
1968-04-09 * | 02:29 | 33.2 | -116.1 | 6.4
1954-03-19 * | 09:54 | 33.3 | -116.2 | 6.2
1952-07-21 * | 11:52 | 35.0 | -119.0 | 7.2
1948-12-04 * | 23:43 | 33.9 | -116.4 | 6.5
1947-04-10 * | 15:58 | 35.0 | -116.6 | 6.2
1946-03-15 * | 13:49 | 35.7 | -118.1 | 6.3
1937-03-25 * | 16:49 | 33.4 | -116.3 | 6.0
1933-03-11 * | 01:54 | 33.6 | -118.0 | 6.3
1925-06-29 * | 14:42 | 34.3 | -119.8 | 6.3
1923-07-23 * | 07:30 | 34.0 | -117.3 | 6.3
1922-03-10 | 11:21 | 35.8 | -120.3 | 6.5
1918-04-21 * | 22:32 | 33.8 | -117.0 | 6.8
1916-11-10 | 09:11 | 35.5 | -116.0 | 6.1
1910-05-15 | 15:47 | 33.7 | -117.4 | 6.0
1908-11-04 | 08:37 | 36.0 | -117.0 | 6.5
1899-12-25 | 12:25 | 33.8 | -117.0 | 6.4
1899-07-22 | 20:32 | 34.3 | -117.5 | 6.5
1892-05-28 | 11:15 | 33.2 | -116.2 | 6.3
1890-02-09 | 12:06 | 33.4 | -116.3 | 6.3
1883-09-05 | 12:30 | 34.2 | -119.9 | 6.0
1858-12-16 | 10:00 | 34.0 | -117.5 | 6.0
1857-01-09 * | 16:24 | 35.7 | -120.3 | 7.6
1855-07-11 | 04:15 | 34.1 | -118.1 | 6.0
1812-12-21 | 19:00 | 34.2 | -119.9 | 7.1
1800-11-22 | 21:30 | 33.0 | -117.3 | 6.5
Table 2: List of earthquake events with $5.3\leq$ M $<6.0$. The symbol * signifies that the event was taken into account by K & K. Date | hr:min |  lat (∘) |  long (∘) | mag
---|---|---|---|---
| (UTC) |  |  |
2018-04-05 | 19:29 | 33.8 | -119.7 | 5.4
2012-08-26 | 20:57 | 33.0 | -115.5 | 5.5
2008-07-29 | 18:42 | 33.9 | -117.8 | 5.9
2005-06-12 | 15:41 | 33.5 | -116.6 | 5.6
2004-09-28 | 17:15 | 35.8 | -120.4 | 5.7
1995-09-20 | 23:27 | 35.8 | -117.6 | 5.8
1995-08-17 | 22:39 | 35.8 | -117.7 | 5.4
1991-06-28 | 14:43 | 34.3 | -118.0 | 5.9
1990-02-28 | 23:43 | 34.1 | -117.7 | 5.5
1988-06-10 | 23:06 | 34.9 | -118.7 | 5.4
1987-10-01 | 14:42 | 34.1 | -118.1 | 5.9
1986-07-08 | 09:20 | 34.0 | -116.6 | 5.6
1981-09-04 | 15:50 | 33.7 | -119.1 | 5.3
1981-04-26 | 12:09 | 33.1 | -115.6 | 5.6
1980-02-25 * | 10:47 | 33.5 | -116.5 | 5.5
1973-02-21 * | 14:45 | 34.1 | -119.0 | 5.9
1970-09-12 * | 14:30 | 34.3 | -117.5 | 5.4
1969-04-28 * | 23:20 | 33.3 | -116.3 | 5.8
1961-01-28 * | 08:12 | 35.8 | -118.0 | 5.3
1955-12-17 * | 06:07 | 33.0 | -115.5 | 5.4
1954-01-12 * | 23:33 | 35.0 | -119.0 | 5.9
1950-07-28 * | 17:50 | 33.1 | -115.6 | 5.4
1949-05-02 * | 11:25 | 34.0 | -115.7 | 5.9
1947-07-24 * | 22:10 | 34.0 | -116.5 | 5.5
1946-07-18 * | 14:27 | 34.5 | -116.0 | 5.6
1946-01-08 * | 18:54 | 33.0 | -115.8 | 5.4
1945-08-15 * | 17:56 | 33.2 | -116.1 | 5.7
1945-04-01 * | 23:43 | 34.0 | -120.0 | 5.4
1944-06-12 * | 11:16 | 34.0 | -116.7 | 5.3
1943-12-22 * | 15:50 | 34.3 | -115.8 | 5.5
1943-08-29 * | 03:45 | 34.3 | -117.0 | 5.5
1941-11-14 * | 08:41 | 33.8 | -118.3 | 5.4
1941-07-01 * | 07:50 | 34.4 | -119.6 | 5.9
1940-05-18 * | 05:03 | 34.1 | -116.3 | 5.4
1938-05-31 * | 08:34 | 33.7 | -117.5 | 5.5
1933-10-02 * | 09:10 | 33.8 | -118.1 | 5.4

To analyse the statistical relationship between lunisolar tidal cycles and
earthquake events, we use the same $\chi^{2}$ algorithm as in the work of K &
K. It is based on the significance of a cyclic Poisson process with respect to
a homogeneous Poisson process evaluated starting from the hypothesis of a
cyclic behaviour with a known phase shift. It enables us to prove the presence
of a clustering of events around a specific phase value with respect to a
specific oscillation with a given period. Of course this method assumes that
the triggering of earthquakes might occur preferentially in some part of a
fundamental cycle characterizing the earth tides. If $N$ is the total number
of events which are being considered, then the $\chi^{2}$ parameter is defined
as :

$\chi^{2}=\frac{2}{N}\left(\sum^{n}_{i=1}\cos{\Phi_{i}}\right)^{2}.$ (1)

Here $\Phi_{i}=2\pi t_{i}/T$ is the phase, where $t_{i}$ is the time of
occurrence of the $i$th. earthquake relative to the time of maximum possible
correlation and $T$ is the period of the oscillation. A zero phase
characterizes an earthquake occurring at a peak of the cycle. A null
hypothesis can be rejected if $\chi^{2}>3.84$. In other words, such a value of
$\chi^{2}$ means that there is a correlation between the lunisolar cycle
considered and the earthquakes at a 95% confidence level. In the same manner,
a value of $\chi^{2}>2.71$ validates a correlation at a 90% confidence level.

We focus our analysis into five lunisolar cycles. First, we studied the
correlation with the average time of interval between rise and fall of the
sun, hereafter named ”semi-diurnal” period. Here the zero phase is taken at 6
a.m and 6 p.m local time and the oscillation period is 12 hrs. Second, we
consider the cycle of the precession of the line of nodes of the Moon orbit
with a 18.613 years period and corresponding by far to the most significant
tidal oscillation at low frequency. Here the zero phase is taken on January
12th., 1932. We also estimated the correlation with other lunisolar cycles
such as the mean fortnightly component (period = 14.77 d.; zero phase =
full/new moon), the mean semi-lunar day (period: 12.42 hrs.; zero phase =
lunar rise/set), and the annual cycle (period = 1 y.; zero phase = December
21th.).

##  3 Phase analysis

Figure 1: Phase distribution of earthquakes occurrence ($\geq 6$ M) with
respect to each lunisolar cycle. Zero phase corresponds to the maximum peak of
the cycle. Here $\bullet$ represents earthquake taken into account by K & K.
Figure 2: Phase distribution of earthquake $5.3\leq M<6$ with respect to each
lunisolar cycle. Zero phase corresponds to the maximum peak of the cycle. Here
$\blacktriangle$ represents the earthquakes taken into account by K & K.

Phase distribution can be used as a rough approach to identify the clustering
in the earthquake events. To calculate the phase $\Phi$, we use the relation
$\Phi_{i}=2\pi t_{i}/T$ and we re-scale the phase to the range (-0.5, 0.5).
Then the phase distribution is classified into two groups: the first one
corresponds to $-0.25\leq\Phi_{i}\leq 0.25$ while the second one is outside
this range.

Figures 1 and 2 show the phase distribution of the earthquake events relative
to each lunisolar cycle considered respectively for the earthquakes with
amplitude M $\geq 6$ and for those with $5.3\leq M<6$. For the precession
cycle, we remark that strong earthquakes are concentrated in the central half
of the diagram with 18 events inside against 12 events outside while weak
earthquakes are more evenly distributed. A similar inequality also appears in
more prominent manner in the case of the semi-diurnal cycle, with 21
earthquakes inside the central half-zone against only 9 outside. Notice that
clustering for both cycles has been recognized by K & K even with a noticeably
smaller number of events taken into account. Indeed they found that from 13
large events recorded with $M\geq 6$, 10 of them were located inside the
central zone for the precession cycle and 9 events inside in the case of semi-
diurnal. Meanwhile, for other cycles, the events are more evenly distributed
both for strong and weak earthquakes.

(a)

(b)

Figure 3: Phase versus phase of two cycles. Here $\bullet$ corresponds to the
earthquake $\geq 6$ M while the earthquake $<6$ is represented by
$\blacktriangle$.

Furthermore, we analyse the possibility of a link between two cyclic tidal
components. The underlying idea is that events occur preferentially (or on the
opposite unfavourably) when they are located simultaneously inside two half-
zones of both cycles. This hypothesis was demonstrated for instance by Klein
(1976) who showed that there is a significant tendency for earthquakes in
selected regions of the world to occur at times at which the tidal stresses
are so aligned as to provide maximum enhancement of the stresses associated
with plate tectonics. For the purpose, we represent a bi-axial plot of the two
phase distributions in which we divide the diagram into four equal quadrants.
As shown in Figure 3, for the case of precession phase versus the semi-diurnal
phase (Figure 3(a)), there is a neat concentration in the upper left zone for
earthquakes with M $\geq 6$. On the contrary, only one earthquake is found in
the upper right zone. Note that such a probability $P$ to find only one event
in one quadrant in the case of an even distribution is:

$P=N\times{\frac{1}{4}}\times\left(\frac{3}{4}\right)^{N-1}$ (2)

where $N$ is the total number of earthquakes. By applying this formula with
$N=32$ in our case we find $P=1.0715\times 10^{-3}\approx 1/933$. This very
small value leads in favour of a bi-modal anti-correlation characterized by a
very high probability that no event occurs when both phases (for solar-day and
precession cycles) lie in the interval [0 - 0.5].

We remark that this trend also appears, but less neatly, when we plot the
phase of the semi-lunar day vs. the phase of the fortnightly component (Figure
3(b)), with only four earthquakes in the upper right zone. Meanwhile, other
bi-modal phase vs. phase configurations do not follow this trend and the
earthquakes are more equivalently distributed.

##  4 $\chi^{2}$ test of correlation

In this section we apply the $\chi^{2}$ statistical test presented in Eq. 1 to
measure the presence of clustering in a quantitative way. For that purpose, we
consider three kinds of data for estimating the $\chi^{2}$ value. First, we
consider all the data starting from 1800 and restricting only for the 32
earthquakes with M $\geq 6$. Second, we consider relatively accurate data
starting from 1933 and we still limit our set to the 15 earthquakes with M
$\geq 6$. Third, we compute the $\chi^{2}$ by considering the weaker
earthquakes from 1933 with magnitude $5.3\leq M<6$. This choice is made
necessary because some data before 1933 is not always accurate and presents a
lack of homogeneity (Kilston and Knopoff, 1983).

Table 3: Value of the $\chi^{2}$ parameter by considering different time range as well as different magnitude threshold. Cycle | Period | 1800 - 2020 | 1933 - 2020 | 1933 - 2020
---|---|---|---|---
|  | $M\geq 6.0$ | $M\geq 6.0$ | $5.3\leq M<6.0$
Precession | 18.613 y. | 1.92 | 3.03 | 0.03
Semi-diurnal | 12 h. | 2.26 | 7.62 | 0.17
Semi-lunar day | 12.421 h. | 0.15 | 1.22 | 0.01
Fortnightly | 14.765 d. | 0.42 | 2.16 | 2.24
Annual | 1 y. | 0.03 | 0.01 | 3.93
No. of events |  | 32 | 15 | 36
Table 4: Comparison of $\chi^{2}$ parameter from this study for the interval [1933-2020] with the one obtained by Kilston and Knopoff (1983) (K & K) for the interval [1933-1980]. Cycle | Period | $M\geq 6.0$ |  | $5.3\leq M<6.0$
---|---|---|---|---
|  | This study | K & K |  | This study | K & K
Precession | 18.613 y. | 3.03 | 5.58 |  | 0.03 | 0.28
Semi-diurnal | 12 h. | 7.62 | 5.20 |  | 0.17 | 0.49
Semi-lunar day | 12.421 h. | 1.22 | 0.35 |  | 0.01 | 0.03
Fortnightly | 14.765 d. | 2.16 | 8.49 |  | 2.24 | 1.23
Annual | 1 y. | 0.01 | 0.05 |  | 3.93 | 0.46
No. of events |  | 15 | 9 |  | 36 | 22

The results of our $\chi^{2}$ tests are displayed in Table 3. As mentioned in
Sect. 2, in the same way as K & K, we accept the null hypothesis for an
homogeneous Poisson process at the 95 % confidence level if $\chi^{2}<3.84$,
and at the 90 % confidence level if $\chi^{2}<2.71$. For the data starting
from 1800, there is no significant correlation between any of the five cycles
considered and the earthquake events. However, if we limit our investigations
to the data from 1933, with homogeneous and accurate values of the magnitude,
we note both a significant correlation of strong earthquakes with precession
and a very big correlation with the semi-diurnal component ($\chi^{2}=7.62$).
On the contrary, the $\chi^{2}$ value is close to zero when considering the
weaker earthquakes dataset both for precession and semi-diurnal components,
which tends to signify a lack of correlation.

Surprisingly, the semi-lunar day component leads to small $\chi^{2}$ values
whatever be the configuration, though the corresponding tidal amplitude is
significantly larger than the semi-diurnal one. Note that the same paradoxical
conclusion was already mentioned by K & K who initially expected that a
correlation should appear with the semi-lunar day cycle. Their argument is
that the lunar tide is considerably stronger than the solar tide, and if a
correlation is found with the semi-diurnal component, an equally strong
correlation should be found with the Moon rise / set periodicity. The reason
of an opposite result might be caused by the chosen period: indeed the lunar
day is not constant, for we adopted a mean value. So the level of correlation
found should strongly depend on the adopted value. This problem is fully
explored with more details in Section 5. The same remarks are available for
the fortnightly cycle whatever be the configuration. Concerning the annual
component, there is no correlation when considering the strong earthquakes
dataset. Meanwhile, taking the earthquakes between 5.3 and 6.0 leads to a
possible correlation ($\chi^{2}=3.93$).

Table 4 shows the comparison of our results with those of K & K whose dataset
did not take into account the earthquakes since 1983. Even if the $\chi^{2}$
values are something different for precession and semi-diurnal components,
they are in both cases above the threshold which leads to the conclusion of
the existence of a correlation for the strong earthquakes. On the opposite the
lack of correlation with the semi-lunar day component is observed in both
studies whatever be the dataset considered.

At last, both studies disagree on the fortnightly and annual components. For
the first one, we find no significant correlation ($\chi^{2}=2.16$) whereas K
& K found a high correlation for strong earthquakes ($\chi^{2}=8.49$). This
could be due to the methodology used : whereas we adopted a mean value of the
period of the fortnightly cycle, K & K referred the occurrence of the
earthquakes to the times of the extrema, rather than to perfectly periodic
occurrence. Concerning the annual cycle, we found a correlation with weak
earthquakes ($\chi^{2}=3.93$) which was not reckoned by K & K. This strongly
raises the question of the variability of the results according to the choice
of the time interval chosen for the dataset.

##  5 Sensitivity test

Following Eq. 1, it is clear that the $\chi^{2}$ value depends on the period
chosen for the corresponding cycle. Indeed, for the lunar day and the
fortnightly cycles, the period between two tidal maxima (or minima) varies in
a significant way, whereas we chose their sole mean value to estimate
$\chi^{2}$ in the previous section. Therefore, here we test the sensitivity of
the $\chi^{2}$ value to the period of the cycle considered by re-estimating
this parameter through the interval between its minimum and maximum values of
the period. We concentrate our study on the two cycles above. Notice that the
lunar month has an average period of about 29.53 days and varies approximately
within the interval $\pm$ 7h. So, the minimum value of the fortnightly period
is 14.59 d. whereas the maximum one is 14.91 d. Meanwhile, the variation of
semi-lunar day is about $\pm$ 9 mn. Hence, the minimum and maximum periods are
respectively 12.27 h. and 12.57 h. Here we restricted the calculations by
using the earthquakes from 1933 and $\geq 6$ magnitude, for the sake of
homogeneity.

As shown in Figure 4(a), the value of $\chi^{2}$ concerning the semi-lunar day
correlation varies between 0 (period = 12.39 hours) and 16 (period = 12.33
hours), thus signifying a high sensitivity with respect to the period
considered. This statement is emphasized by the presence of numerous peaks in
the interval.

To illustrate the high level of correlation found for period 12.33 h. we show
in Figure 4(b) the corresponding phase distribution of the earthquakes. For
the 15 earthquakes with M $\geq 6$, the phase distribution is much more
centralized in the middle half of the phase. Only one earthquake is found
outside. By adopting the same formula as Eq. 2 with N=15, we obtain a
probability of P = 6.68 % that such a distribution is obtained by considering
a null hypothesis. However, we remark that such a clustering besides being
still apparent, is less obvious for earthquakes with M $<6$.

(a) $\chi^{2}$ semi-lunar day cycle

(b) ”12.33 hours” phase

(c) $\chi^{2}$ fortnightly cycle

(d) ”14.86 days” phase

Figure 4: Left: $\chi^{2}$ estimation between the minimum and maximum values
of the semi-lunar day period (a) and fortnightly (c). Red dashed line
represents the nominal value that we used in Section 4. Right: phase
distribution with respect to the phase of the semi-lunar day with period 12.33
hours (b) and the phase of fortnightly with period 14.86 days (d). Here
$\bullet$ represents the earthquake with M $\geq 6$ and $\blacktriangle$
represent the earthquake with M $<6$.

Following a similar procedure, Fig. 4(c) depicts the variation of $\chi^{2}$
in the range of possibilities of the value of the fortnightly period. Here
$\chi^{2}$ varies from almost 0 (for a 14.81 d. period) to a rather limited
value of 6.5 (for 14.86 d.). Figure 4(d) shows the phase distribution with the
period = 14.86 d. corresponding to the maximum. In a similar way to Fig. 4(b),
the earthquake phase is located in the middle half of the figure for the
earthquakes with M $\geq 6$. Accordingly we suggest being more cautious for
interpreting our results related to $\chi^{2}$ from semi-lunar day and
fortnightly periodicities. A more rigorous method needs to be applied to study
the correlation between these two cycles and earthquakes events, in particular
by computing the phase by interpolation with respect to two successive extrema
of the cycle.

Furthermore, it is worthwhile to analyze the phase vs. phase relation between
the maximum $\chi^{2}$ period of semi-lunar day (period = 12.33 hr.) and of
the fortnightly period (14.86 d). As shown in Fig. 5, the earthquakes with M
$\geq 6$ are almost evenly distributed among four different quadrants.
However, they are significantly concentrated in the centre of the plot. It
likely signifies that the probability that an earthquake occurs becomes higher
when the two events are at their peak (zero phase). Besides these two cycles,
there is no meaning to carry out such a $\chi^{2}$ sensitivity test for
precession, semi-diurnal and annual cycles, for they do not exhibit
significant period variability.

Figure 5: Earthquake phase distribution in fortnightly with period = 14.86 d.
vs. the phase distribution in semi-lunar day with period 12.33 h. Earthquakes
with M $\geq 6$ are represented in $\bullet$,those with M $<6$ M are in
$\blacktriangle$.

##  6 Conclusion

In this study, we re-conducted a work carried out nearly four decades ago by K
& K to investigate the possibilities of correlation between some lunisolar
cycles and earthquake events that occurred in the Southern California region.
We beneficiate of a considerably larger and more accurate data than in this
last study. First, we analysed the phase distribution of the earthquakes with
respect to the leading lunisolar cyclic tidal components. The corresponding
periods are the semi-diurnal one (12 h.), the semi-lunar day (12.42 d) the
fortnightly period(14.77 d), the year (365.25 d) and the period of the
precession of the nodes of the orbit of the Moon (18.613 y.). We remarked that
apparent clustering of earthquakes with M $\geq 6$ occurs near the zero phase
in the diagrams related to the precession and semi-diurnal cycles, thus
confirming a statement already established by K & K. Moreover, we analysed the
possible link between these two cycles. We found that there is a high
probability that no event occurs when both phases (for semi-diurnal and
precession cycles) lie in the interval of [0, 0.5], with only a 0.11 %
probability of the studied distribution to happen when adopting the hypothesis
of a lack of correlation.

In order to confirm the tendencies found above, we applied the same $\chi^{2}$
statistical test as that established by K & K to validate the presence of
clustering in a more rigorous way. We found that when limiting our
investigations to events between 1933 and 2020 for which the quality of data
is quite homogeneous and particularly accurate, strong earthquakes with M
$\geq 6$ present a high level of correlation with semi-diurnal and precession
phase, which clearly confirms K & K conclusions. In addition, we suspect also
the presence of a significant correlation between weaker earthquakes ($5.3\leq
M<6.0$) and the annual cycle, which was not pointed out by K & K. At last
there is apparently no correlation of earthquakes with the fortnightly and
semi-lunar day tidal cycles, despite their strong respective tidal amplitudes.
To measure the robustness of our result, we re-estimated $\chi^{2}$ by varying
the period of the fortnightly and semi-lunar day cycles between their minimum
and maximum values. Thus we found that the $\chi^{2}$ results for both cycles
are not robust. In other words, the value of $\chi^{2}$ can assert a total
absence or a high level of correlation according to the value of the period
selected. Therefore a more rigorous statistical method is needed to study
these two cycles, which should consider phase shifts due to the variability of
the periods considered. These phase shifts, if not taken into account, would
lead to misidentification of the times of peaks of tidal components, as
already explained by K & K. Such a sensitivity test is not necessary for the
annual, diurnal (12h.), and precession components, which have a well
established and not varying period.

An important output from the present study is that we confirm the main results
found by K & K a little less than four decades ago concerning the interaction
of luni-solar extensional Earth tides in a part of the Earth where strike-slip
faulting occurs with a well-defined geographical orientation. It should be
interesting to extrapolate our positive results concerning a tides-earthquakes
correlation for the semi-diurnal and precession cycles to anticipate the time
of the next episode of high potential for the occurence of large earthquakes
in southern California.

## References

  * Cadicheanu, Van Ruymbeke, and Zhu (2007) Cadicheanu, N., Van Ruymbeke, M., Zhu, P.: 2007, Tidal triggering evidence of intermediate depth earthquakes in the Vrancea zone (Romania). Nat Hazard Earth Syst Sci 7, 733 – 740.
  * Center (2013) Center, S.C.E.D.: 2013, Southern California Earthquake Center. Caltech. Dataset.
  * Chen et al. (2012) Chen, H.J., Chen, C.C., Tseng, C.Y., Wang, J.H.: 2012, Effect of tidal triggering on seismicity in Taiwan revealed by the empirical mode decomposition method. Natural Hazards and Earth System Sciences 12(7), 2193\.
  * Kilston and Knopoff (1983) Kilston, S., Knopoff, L.: 1983, Lunar–solar periodicities of large earthquakes in southern California. Nature 304(5921), 21 – 25.
  * Klein (1976) Klein, F.W.: 1976, Earthquake swarms and the semidiurnal solid earth tide. Geophysical Journal International 45(2), 245 – 295.
  * Knopoff (1964) Knopoff, L.: 1964, The statistics of earthquakes in Southern California. Bulletin of the Seismological Society of America 54(6A), 1871 – 1873.
  * Métivier et al. (2009) Métivier, L., de Viron, O., Conrad, C.P., Renault, S., Diament, M., Patau, G.: 2009, Evidence of earthquake triggering by the solid earth tides. Earth and Planetary Science Letters 278(3-4), 370 – 375.
  * Moncayo, Zuluaga, and Monsalve (2019) Moncayo, G.A., Zuluaga, J.I., Monsalve, G.: 2019, Correlation between tides and seismicity in Northwestern South America: The case of Colombia. Journal of South American Earth Sciences 89, 227 – 245.
  * Nikolaev (1995) Nikolaev, V.: 1995, On the earthquake triggering by the Earth’s tide and underground nuclear explosions. In: Earthquakes Induced by Underground Nuclear Explosions, Springer, ???, 255 – 265.
  * Palumbo (1986) Palumbo, A.: 1986, Lunar and solar tidal components in the occurrence of earthquakes in Italy. Geophysical Journal International 84(1), 93 – 99.
  * Ryall, VanWormer, and Jones (1968) Ryall, A., VanWormer, J., Jones, A.E.: 1968, Triggering of microearthquakes by earth tides, and other features of the Truckee, California, earthquake sequence of September, 1966. Bulletin of the Seismological Society of America 58(1), 215 – 248.
  * Saltykov, Ivanov, and Kugaenko (2004) Saltykov, V., Ivanov, V., Kugaenko, Y.A.: 2004, The action of Earth tides on seismicity prior to the Mw 7.0 November 13, 1993, Kamchatka earthquake. Fizika Zemli 40(7), 563 – 572.
  * Shirokov (1983) Shirokov, V.: 1983, The influence of the 19-year tidal cycle on large-scale eruptions and earthquakes in Kamchatka, and their long-term prediction. The Great Tolbachik Fissure Eruption, 232 – 241.
  * Shlien (1972) Shlien, S.: 1972, Earthquake-tide correlation. Geophysical Journal International 28(1), 27 – 34.
  * Shudde and Barr (1977) Shudde, R., Barr, D.: 1977, An analysis of earthquake frequency data. Bulletin of the Seismological Society of America 67(5), 1379 – 1386.
  * Souchay and Stavinschi (1997) Souchay, J., Stavinschi, M.: 1997, Study of the correlations between long-periodic terrestrial tides and occurrence of earthquakes in the Vrancea site. Earth, Moon, and Planets 77(2), 105 – 124.
  * Stavinschi and Souchay (2003) Stavinschi, M., Souchay, J.: 2003, Some correlations between earthquakes and earth tides. Acta Geodaetica et Geophysica Hungarica 38(1), 77 – 92.
  * Stover and Coffman (1993) Stover, C.W., Coffman, J.L.: 1993, Seismicity of the United States, 1568-1989 (revised), US Government Printing Office, ???.
  * Tamrazyan (1974) Tamrazyan, G.: 1974, Possible cosmic influences on the 1966 Tashkent earthquake and its largest aftershocks. Geophysical Journal International 38(2), 423 – 429.
  * Tanaka, Ohtake, and Sato (2002) Tanaka, S., Ohtake, M., Sato, H.: 2002, Evidence for tidal triggering of earthquakes as revealed from statistical analysis of global data. Journal of Geophysical Research: Solid Earth 107(B10), ESE1.
  * Vergos, Arabelos, and Contadakis (2015) Vergos, G., Arabelos, D., Contadakis, M.: 2015, Evidence for Tidal triggering on the earthquakes of the Hellenic Arc, Greece. Physics and Chemistry of the Earth, Parts A/B/C 85, 210 – 215.
  * Young et al. (1978) Young, D., Zurn, W., et al.: 1978, Tidal triggering of earthquakes in the Swabian Jura? Journal of Geophysics 45(1), 171 – 182.


# Encoding Domain Knowledge in Multi-view Latent Variable Models:
A Bayesian Approach with Structured Sparsity

Arber Qoku  German Cancer Consortium (DKTK) and German Cancer Research Center
(DKFZ)  Goethe University Frankfurt, Germany  Florian Buettner  German Cancer
Consortium (DKTK) and German Cancer Research Center (DKFZ)  Goethe University
Frankfurt, Germany  Frankfurt Cancer Insititute, Germany
{arber.qoku<EMAIL_ADDRESS>
###### Abstract

Many real-world systems are described not only by data from a single source
but via multiple data views. For example, in genomic medicine, a patient can
be described by data from different molecular layers. This raises the need for
multi-view models that are able to disentangle variation within and across
data views in an interpretable manner. Latent variable models with structured
sparsity are a commonly used tool to address this modeling task but
interpretability is cumbersome since it requires a direct inspection and
interpretation of each factor via a specialized domain expert. Here, we
propose MuVI, a novel approach for domain-informed multi-view latent variable
models, facilitating the analysis of multi-view data in an inherently
explainable manner. We demonstrate that our model (i) is able to integrate
noisy domain expertise in form of feature sets, (ii) is robust to noise in the
encoded domain knowledge, (iii) results in identifiable factors and (iv) is
able to infer interpretable and biologically meaningful axes of variation in a
real-world multi-view dataset of cancer patients.

##  1 Introduction

Figure 1: Performance of our model (MuVI) against competitive models. Only our
model achieves a higher $F_{1}$ score for recovering true active features
while maintaining a low reconstruction error.

In many real-world applications, complex systems are characterized via
multiple data views. That is, observations are represented by multiple groups
of distinct features. These groups of features often represent different and
complementary input sources that are required for a comprehensive
characterization of a sample. For example, in genomic medicine, a patient can
be described by quantifying different molecular layers such as the proteome,
the microbiome and the transcriptome.
Latent variable models are powerful statistical tools that uncover the axes of
variation between samples and data-views, by inferring unobserved hidden
states from the observable high-dimensional data. To disentangle the sources
of heterogeneity driving intra-view and inter-view variation in a meaningful
manner, it is key that latent variable models are expressive and
interpretable. Modeling approaches with the highest expressive power are based
on variational autoencoders[Ainsworth et al., 2018] or Gaussian Process Latent
Variable models [Damianou et al., 2012]. However, the non-linearities
introduced to ensure a high expressive power hinder interpretability to the
extent that linear approaches remain the dominant tool to analyze complex
multi-view data in many domains[Argelaguet et al., 2018]. A recently proposed
AE-based approach attempts to balance expressive power and interpretability by
combining a non-linear encoder with a linear decoder[Svensson et al., 2020].
The most commonly used modeling framework for yielding an interpretable
decomposition of multi-view data are factor analysis models with structured
sparsity [Zhao et al., 2016, Klami et al., 2014]. These models result in
sparse factor loadings, allowing for a direct inspection and interpretation of
each factor from a domain expert.
However, this manual annotation of factors is cumbersome, time consuming and
requires highly-specialized expertise.
To automate this annotation process, we propose to leverage the partial
knowledge on the structure of the latent space that is readily available in a
structured manner in many domains. For example, in genomic medicine, disease
states are often characterized by sets of features (or pathways), comprising
genes that are known to act in a coordinated manner. These data, collected
over many years of scientific progress, are curated in dedicated pathway
databases.
However, integrating this domain knowledge in a principled manner as sparsity
prior in multi-view models is challenging: domain knowledge is noisy (e.g.
pathways contain many false positive and false negative annotations) and
domain knowledge is incomplete and often only available for a subset of views
(e.g. in genomic medicine pathway information is only readily available for
the transcriptome and the genome).
In this work we propose a novel approach for incorporating domain knowledge
via structured sparse priors for analysing multi-view data in an inherently
explainable manner. Our domain-informed priors render the learnt latent
variables directly interpretable without any further interaction with domain
expert, by tagging each latent variable with its corresponding feature set (or
pathway in genomic medicine).
The contributions of our paper are as follows:

  * •

We develop a multi-view latent variable model with domain-informed structured
sparsity (MuVI) that integrates noisy domain expertise in terms of feature
sets.

  * •

We demonstrate that our model utilizes prior information efficiently by
recovering correct signals from noisy feature sets and is robust to poorly
specified priors.

  * •

We show that feature-set priors result in identifiable factors, and promote
partial identifiability in the uninformed views.

  * •

We demonstrate the practical utility of our model on a large multi-view
dataset of cancer patients by inferring interpretable and biologically
meaningful axes of variation.

##  2 Related Work

Factor analysis (FA) is a fundamental approach for estimating and
understanding the correlation structure among observed variables[Thurstone,
1931], which has inspired the development of numerous latent variable models.
Due to its simplicity, however, standard FA is unable to model observations
from multiple sources. Extensions, such as the canonical correlation analysis
(CCA)[Hotelling, 1992, Klami et al., 2013], or the group factor analysis
(GFA)[Klami et al., 2014] model paired observations simultaneously by learning
linear dependencies underlying two or more data sources or views. A central
assumption of the GFA is that the multi-view observations are a manifestation
of a lower dimensional common latent space corrupted by Gaussian noise. The
main components that describes this relationship are the factor-to-feature
linear mappings or the so-called factor loadings. Factor loadings encode the
structure of each factor, and play an important role in the interpretation of
the model. Hence, factor loadings are suitable for introducing statistical
assumptions on their underlying structure such as sparsity. For instance, GFA
quantifies the association between view and factor by extending the automatic
relevance determination (ARD)[MacKay et al., 1994]. In non-Bayesian
approaches, sparsity is commonly handled by introducing additional terms of an
optimization objective posed on the model parameters. A canonical example is
the the L1 penalty, also known as the lasso[Tibshirani, 1996]. On the other
hand, Bayesian approaches achieve sparse solutions via particular sparsity
inducing priors. For instance, the double exponential, or Laplace prior is the
Bayesian counterpart of the lasso[Park and Casella, 2008]. Another popular
choice is the discrete spike-and-slab prior[Mitchell and Beauchamp, 1988], a
mixture of a Dirac delta distribution centered at zero for pruning irrelevant
signals, and a normal distribution for modeling larger signals. Recently, the
spike-and-slab lasso (SSL)[Ročková and George, 2018] has emerged as a
combination of two Laplace distributions that are parameterized to emulate
both the spike and the slab component. Other flavours of continuous shrinkage
priors[Polson and Scott, 2010] include the horseshoe prior[Carvalho et al.,
2009, 2010], which offers a continuous relaxation of the spike-and-slab
approach, thereby preserving sparse properties while providing computational
benefits such as being differentiable.

Several Bayesian approaches have successfully utilized the synergy between
latent variable models and sparsity inducing priors[Bernardo et al., 2003,
Engelhardt and Stephens, 2010, Lan et al., 2014]. In the multi-view setting,
[Zhao et al., 2016] propose a hierarchical Bayesian approach for group factor
analysis with structured sparsity, facilitated by a cascading three parameter
beta prior[Armagan et al., 2011]. The consequences of this addition are
twofold, the column-wise sparsity supports the inference of associations
between views, while element-wise sparsity encourages feature selection within
each individual factor. Multi-omics factor analysis (MOFA)[Argelaguet et al.,
2018] assumes a similar structured sparsity in two levels, which is achieved
by combining an ARD prior with a spike-and-slab prior. Beyond linear models,
there is an array of non-linear modeling approaches that generalize the ideas
of linear LVMs and sparse priors to more complex data settings. These include
approaches based on variaitonal autoencoder [Ainsworth et al., 2018] or
Gaussian Process Latent Variable Models [Damianou et al., 2012]. Such
approaches, however, are difficult to interpret in practice and have therefore
not received widespread application. For example, in the latter GP-based
approach, loading matrices are marginalized out, and cannot be used to
identify any physical or biological processes captured by individual factors.
To maintain the expressive power of autoencoders and still facilitate the
interpretability of linear models, [Svensson et al., 2020] have recently
proposed to combine a non-linear encoder with a linear decoder in an
autoencoder framework. While this approach was proposed for single data views,
extension to the multi-view setting is straight-forward.
Applications of shrinkage priors are not limited to feature selection, and
have been successfully applied in other supervised tasks such as model
selection in Bayesian neural networks[Ghosh et al., 2019, Louizos et al.,
2017].

Typically, before training sparse latent variable models, the prior
assumptions are set to be relatively vague and uninformative, allowing only
the observations to mold the posterior. An alternative approach is introducing
an informed structured sparsity into the model based on external expertise in
the domain of interest. The factorial single-cell latent variable model
[Buettner et al., 2017] attempts to bridge this gap by explicitly integrating,
and jointly modeling a collection of domain relevant feature sets such as gene
set annotations. The model incorporates and extends the spike-and-slab to be
informed by the presence or absence of a feature, and infer pre-labeled latent
axes of variation. However, this model is only able to model a single data
view and the proposed inference scheme is cumbersome and does not scale to
large sample sizes. An orthogonal approach for identifying latent variables
corresponding to pathways is the single sample gene set enrichment analysis
(ssGSEA) [Barbie et al., 2009], where an absolute enrichment of a pathway in
each sample is computed. However, this approach does not account for noisy
annotations, treats all pathways as independent and can only be applied to a
single data view.
A multi-view approach with structured sparsity that admits prior information
in terms of feature sets, and is able to handle partial and potentially noisy
priors has not yet been explored.

##  3 Domain-informed multi-view modeling

###  3.1 Background and notation

Let $\mathbf{y}_{i}\in\mathbb{R}^{D}$ denote a $D$-dimensional observation for
$i\in\left\\{1,\dots,N\right\\}$, and $G_{m}\subseteq\\{1,\dots,D\\}$ describe
a mutually disjoint grouping of the features into $M$ data views, where
$G_{p}\cap G_{q}=\emptyset$ for $p\neq q\in\\{1,\dots,M\\}$. For simplicity,
we assume a permutation of the features $D$ such that the first $D_{1}$
features belong to $G_{1}$, the second $D_{2}$ features to $G_{2}$ and so on.
As a result, we may rewrite the collection of observations as a matrix
$\mathbf{Y}\in\mathbb{R}^{N\times D}$ comprising $M$ coupled views
$\left[\mathbf{Y}^{(1)},\mathbf{Y}^{(2)},\dots,\mathbf{Y}^{(M)}\right]$, where
$\mathbf{Y}\in\mathbb{R}^{N\times D_{m}}$. The main goal is to then represent
each observation $\mathbf{y}_{i}$ in terms of a low-dimensional set of latent
factors $\mathbf{x}_{i}\in\mathbb{R}^{K},K\ll D$. The relationship between the
observations and the latent factors is described by a set of view-specific
factor loadings $\mathbf{W}^{(m)}\in\mathbb{R}^{D_{m}\times K}$. Then, the
recipe for the group factor analysis framework posits the following generative
process to the observed data,

$\mathbf{y}_{i}^{(m)}\sim\mathcal{N}\left(\mathbf{W}^{(m)}\mathbf{x}_{i},\mathbf{\Psi}^{(m)}\right)$
(1)

where each latent variable typically follows an isotropic standard normal
distribution

$\mathbf{x}_{i}\sim\mathcal{N}\left(\mathbf{0},\mathbf{I}\right).$ (2)

The residuals are denoted by
$\mathbf{\Psi}^{(m)}=diag(\mathbf{\sigma}^{2(m)})$, a diagonal matrix storing
the marginal variances $\sigma_{j}^{2(m)}$ of each variable $j$ in view $m$.
Due to conjugacy properties, setting

$\sigma_{j}^{2(m)}\sim\Gamma^{-1}(\alpha_{\sigma},\beta_{\sigma})$ (3)

is a common choice, where $\Gamma^{-1}(\alpha,\beta)$ describes the inverse-
Gamma with shape and scale parameters $\alpha$ and $\beta$. Finally, an
important component of GFA is the collection of the factor loadings which
applies a linear projection of the latent variable $\mathbf{x}_{i}$ to
$\mathbf{y}_{i}^{(m)}$. In order to facilitate the interpretability of the
factor-to-feature mapping, [Klami et al., 2014] propose a structured sparsity
for the columns of $\mathbf{W}^{(m)}$, such that each factor falls into two
distinct categories. A non-zero factor loading vector $\mathbf{w}_{k}^{(m)}$
indicates an active factor $k$, while $\mathbf{w}_{k}^{(m)}=\mathbf{0}$ an
inactive factor $k$ in view $m$. Hence, a factor is either shared across an
arbitrary subset of views, or private to a specific view.

###  3.2 MUVI

We follow a similar approach, inspired by the success of the horseshoe
prior[Carvalho et al., 2009, 2010], and introduce a _view-factor-local_
shrinkage prior on the factor loadings to enable both column-wise and element-
wise shrinkage:

$w_{j,k}^{(m)}\sim\mathcal{N}\left(0,\left(\tau^{(m)}\delta_{k}^{(m)}\lambda_{j,k}^{(m)}\right)^{2}\right),$
(4)

where each scale in the hierarchy follows a positive Cauchy distribution,

$\displaystyle\tau^{(m)}$ $\displaystyle\sim\mathcal{C}^{+}(0,1)$ (5)
$\displaystyle\delta_{k}^{(m)}$ $\displaystyle\sim\mathcal{C}^{+}(0,1)$ (6)
$\displaystyle\lambda_{j,k}^{(m)}$ $\displaystyle\sim\mathcal{C}^{+}(0,1).$
(7)

Each level in the hierarchy contributes to the overall structured sparsity of
the factor loadings. In particular, $\delta_{k}^{(m)}$ serves as an automatic
relevance determination (ARD) mechanism[MacKay et al., 1994] for factor $k$ in
view $m$, effectively decoupling this factor from the rest of the views. At
the same time, $\lambda_{j,k}^{(m)}$ acts as a regulator on each individual
loading, encouraging the model to seek simpler solutions that describe each
factor in terms of fewer features. However, due to the heavy tails of the
Cauchy distribution, weakly identified loadings under the horseshoe prior can
easily escape the regularization. To counteract this behaviour, the
regularized horseshoe guarantees a non-zero penalty even for large
weights[Piironen and Vehtari, 2017]. Concurrently, the regularized
parameterization helps integrate prior information about the structure of the
latent factors into the model. We update Equation 4 as follows. To simplify
notation, we sometimes drop the view-specific superscript $m$, which can also
be implicitly encoded in $j$. Let $\gamma_{j,k}=\tau\delta_{k}\lambda_{j,k}$,
then

$w_{j,k}\sim\mathcal{N}\left(0,\frac{\left(c_{j,k}\gamma_{j,k}\right)^{2}}{c_{j,k}^{2}+\gamma_{j,k}^{2}}\right)$
(8)

Note that in Equation 8 the $c_{j,k}^{2}$ parameter relates to a specific
factor loading rather than being a global parameter as is typically the case.
For a relatively large $c_{j,k}^{2}\gg\gamma_{j,k}^{2}$, $w_{j,k}$ is scaled
by a factor of nearly $\gamma_{j,k}$, rendering the effect of the additional
penalty insignificant, and reinstating the original horseshoe prior. On the
other hand, when $c_{j,k}^{2}\ll\gamma_{j,k}^{2}$, $w_{j,k}$ is scaled by a
factor of nearly $c_{j,k}$, thus setting an upper bound on corresponding
factor loadings. Equivalently, the regularized horseshoe can be seen as a
continuous alternative to the discrete spike-and-slab prior[Mitchell and
Beauchamp, 1988] with a finite slab width. The inverse-Gamma is a suitable
choice as a prior

$c_{j,k}^{2}\sim\Gamma^{-1}(\alpha_{c},\beta_{c}).$ (9)

###  3.3 Integrating Prior Knowledge of Feature Sets

Next, we attempt to integrate prior knowledge in terms of potentially noisy
feature sets into our MuVI. A feature set $\mathbf{I}^{D}$ is a collection of
binary variables $I_{j}\in\left\\{0,1\right\\}$,
$j\in\left\\{1,\dots,D\right\\}$, where $I_{j}=1$ indicates the presence of
feature $j$, and $I_{j}=0$ its absence. Assume we have substantial knowledge
about the underlying structure of the factor loadings in terms of feature
sets. That is, for every latent dimension $k$ and a set of features $D$ we are
given a corresponding feature set $\mathbf{I}_{k}^{D}$, such that $w_{j,k}>0$
if $I_{j,k}=1$, and $w_{j,k}=0$ otherwise. Then, the solution to a sparse
group factor analysis is trivial. Each latent factor can be seen as a factor-
to-feature mapping of the subset of active features determined by the feature
set. In terms of Equation 8, this can also be achieved by setting the
corresponding slab width $c_{j,k}\approx 0$, applying an infinitely large
weight decay penalty to $w_{j,k}$. In practice, we rarely have access to such
pristine ground truth. However, in some cases, we may exploit existing domain
knowledge to develop a prior belief about a plausible structure of the factor
loadings by accommodating noisy feature sets into our modeling approach. Let
$\tilde{\mathbf{I}}^{D}$ be a noisy version of a feature set $\mathbf{I}^{D}$,
where a subset of the binary variables
$I_{Q},Q\subset\left\\{1,\dots,D\right\\}$ has been flipped to generate
$\tilde{\mathbf{I}}^{D}$, inserting a non-zero fraction of false positives and
false negatives. In order to successfully integrate noisy feature sets into
MuVI, we need to relax the regularization penalty introduced by a pre-defined
$c_{j,k}$, allowing the adaptation of incorrect signals, given sufficient
evidence from the data. Subsequently, we scale $c_{j,k}$ by an auxiliary
constant $0<\alpha_{j,k}\leq 1$. A value of $\alpha_{j,k}=1.0$ poses no prior
penalty to $w_{j,k}$, while $\alpha_{j,k}<1.0$ leads to a smaller slab width a
priori. Consequently, a smaller $\alpha_{j,k}$ translates to a stronger prior
belief. Experiments show that values around $0.01\leq\alpha_{j,k}\leq 0.05$
for absent features in a prior collection of feature sets provide the best
results on average.

###  3.4 Inference

Finally, the joint model can be written as

$\displaystyle p\left(\mathbf{Y},\mathbf{\Theta}\right)$
$\displaystyle=p\left(\mathbf{Y},\mathbf{X},\mathbf{W},\mathbf{\Lambda},\mathbf{\Delta},\mathbf{\tau},\mathbf{C},\mathbf{\Psi}\right)$
$\displaystyle\times
p\left(\mathbf{Y}\mid\mathbf{X},\mathbf{W},\mathbf{\Psi}\right)p\left(\mathbf{\Psi}\right)$
$\displaystyle\times
p\left(\mathbf{X}\right)p\left(\mathbf{W}\mid\mathbf{\Lambda},\mathbf{\Delta},\mathbf{\tau},\mathbf{C}\right)$
$\displaystyle\times
p\left(\mathbf{\Lambda}\right)p\left(\mathbf{\Delta}\right)p\left(\mathbf{\tau}\right)p\left(\mathbf{C}\right),$
(10)

where
$\mathbf{\Theta}=\mathbf{Y},\mathbf{X},\mathbf{W},\mathbf{\Lambda},\mathbf{\Delta},\mathbf{\tau},\mathbf{C},\mathbf{\Psi}$,
and

$\displaystyle\mathbf{Y}$
$\displaystyle\\!=\\!\left\\{y_{i,j}^{(m)}\right\\},\mathbf{\Psi}\\!=\\!\left\\{\sigma_{j}^{(m)}\right\\},\mathbf{X}\\!=\\!\left\\{x_{i,k}\right\\},\mathbf{W}\\!=\\!\left\\{w_{j,k}^{(m)}\right\\},$
$\displaystyle\mathbf{\Lambda}$
$\displaystyle\\!=\\!\left\\{\lambda_{j,k}^{(m)}\right\\},\mathbf{\Delta}\\!=\\!\left\\{\delta_{k}^{(m)}\right\\},\mathbf{\tau}\\!=\\!\left\\{\tau^{(m)}\right\\},\mathbf{C}\\!=\\!\left\\{c_{j,k}^{(m)}\right\\}.$

Inference is performed by introducing a fully factorized family of
parameterized distributions
$q_{\phi}\left(\mathbf{\Theta}\right)=\prod_{\theta\in\mathbf{\Theta}}q_{\phi}\left(\theta\right)$
to approximate the intractable posterior
$p\left(\mathbf{\Theta}\mid\mathbf{Y}\right)$. We maximize the evidence lower
bound with respect to the variational parameters $\phi$, which, in return
reduces the gap between the true and the approximate posterior in terms of the
KL divergence. The family of normal distributions is a natural choice for
approximating $\mathbf{X}$ and $\mathbf{W}$, whereas, for the rest of the
parameters we assert the log-Normal distribution to ensure positive
samples[Ghosh et al., 2018]. The resulting optimization objective

$\mathcal{L}(\phi)=\mathbb{E}_{q_{\phi}}\left[\log
p(\mathbf{Y},\mathbf{\Theta})-\log
q_{\phi}\left(\mathbf{\Theta}\right)\right]$ (11)

is amenable to black box stochastic variational inference[Ranganath et al.,
2014, Hoffman et al., 2013], which involves sampling from the variational
distribution, and taking unbiased but noisy Monte Carlo estimates of the
gradient $\nabla_{\phi}\mathcal{L}$. In addition, the reparameterization
trick[Kingma and Welling, 2013] further stabilizes the optimization procedure,
in which the random variables are expressed as a combination of deterministic
variables and external random noise, thereby greatly reducing the variance of
the MC estimates. We provide an implementations for optimization using
Pyro[Bingham et al., 2019] on GitHub111https://github.com/MLO-lab/MuVI.

##  4 Experiments

We first evaluate different aspects of MuVI empirically on a wide range of
simulated settings. Ideally, the model performs well in the following general
tasks:

  * •

Learning a meaningful latent representation with low reconstruction loss

  * •

Utilizing prior information efficiently by recovering correct signals from
noisy feature sets

  * •

Being robust against poorly specified or entirely incorrect prior information

  * •

Promoting an implicit flow of the prior information from the informed to the
uninformed views via the shared factors

We compare our model against three baselines that admit observations across
multiple views, and learn a sparse representation of the common latent space:
the Bayesian group factor analysis with structured sparsity (BASS)[Zhao et
al., 2016], multi-omics factor analysis (MOFA)[Argelaguet et al., 2018] and an
interpretable multi-view autoencoder (AE), as a naive multi-view extension to
the interpretable autoencoder[Svensson et al., 2020].
Finally, we demonstrate the practical utility of MuVI on a real-world dataset
from The Cancer Genome Atlas (TCGA)[Tomczak et al., 2015].

###  4.1 Synthetic Experiments

####  4.1.1 Data generation

We compile the synthetic data by generating $N=200$ samples across four views,
each comprising $D_{1}=D_{2}=D_{3}=D_{4}=400$ features. The latent space
consists of $K=15$ factors, that are linearly transformed by a set of sparse
factor loadings. Each weight is sampled independently from a standard normal
distribution, where loadings with an absolute value of less than $0.1$ are set
to zero, to emphasize the gap between active and inactive signals. In
addition, we randomly set $85\%$-$95\%$ of the loadings to zero. The
relationship between the data views is explained by the structured sparsity of
the loadings. In particular, we generate all combinations of four binary
variables, resulting in $15$ distinct relationship configurations for each
factor as depicted in the top heatmap of Figure 2. A dark entry means the
factor does not contribute to the corresponding view, and decouples it from
the rest of the views. For instance, factor with index 0 is fully shared
across all views, whereas factor 7, 11, 13 and 14 explain variability
pertaining to each view individually. The rest of the factors exhibit a
mixture of partially shared configurations.

Figure 2: Depiction of cross-view relationships in the synthetic data. True
relationship (top) versus inferred relationship based on the factor scales
learned by MuVI (bottom).

We extract the true factor loadings mask, indicating active and inactive
features for each latent dimension, and generate potentially noisy feature
sets which serve as prior information during training. We perturb the feature
sets by swapping $10\%,20\%,50\%,90\%$ and $100\%$ of the true positives with
true negatives. The severity of the noise translates to a poorly specified
prior belief, where $100\%$ noise means that the prior belief is entirely
incorrect, up to the number of the expected active features. Also, we
introduce additional false positives for factors that are inactive, as to
further disguise any true underlying relationship within and between views.

####  4.1.2 Evaluation

We evaluate all models on five randomly generated datasets as described above,
and report average metric scores across all views. We assess the quality of
the representation in terms of RMSE between the true observations and the
reconstructed views from the common latent space. We quantify the concordance
between the inferred factor loadings and the true underlying structure by
reporting the precision, recall and the $F_{1}$ score. As a prerequisite to
calculating the binary scores we extract a binary representation of active and
inactive features based on a threshold. For our model and MOFA, the threshold
of $0.1$ matches the true cutoff between active and inactive signals when
generating the sparse loadings. For BASS and AE, the optimal thresholds were
$0.05$ and $0.01$, respectively. Moreover, since the factor analysis suffers
from non-identifiability[Shapiro, 1985], a matching of the latent factors and
their corresponding factor loadings is necessary. However, we highlight that
this issue is less severe for the latent structure learned by our model, and
later show that the structure reinforced by the prior information is, under
certain conditions, identifiable.

| RMSE | Precision | Recall | $F_{1}$
---|---|---|---|---
BASS | $0.543\pm 0.02$ | $\mathbf{0.944\pm 0.05}$ | $0.899\pm 0.03$ | $0.920\pm 0.04$
MOFA | $0.331\pm 0.00$ | $0.541\pm 0.14$ | $0.914\pm 0.02$ | $0.672\pm 0.11$
AE | $0.325\pm 0.01$ | $0.692\pm 0.06$ | $0.964\pm 0.01$ | $0.805\pm 0.04$
MuVI$-$ (uninformed) | $\mathbf{0.323\pm 0.00}$ | $0.933\pm 0.05$ | $\mathbf{0.995\pm 0.00}$ | $\mathbf{0.962\pm 0.03}$
MuVI$0.5$ (1 inf. view) | $\mathbf{0.322\pm 0.00}$ | $0.966\pm 0.01$ | $0.996\pm 0.00$ | $0.981\pm 0.00$
MuVI$0.1$ (1 inf. view) | $\mathbf{0.322\pm 0.00}$ | $0.965\pm 0.01$ | $0.997\pm 0.00$ | $0.982\pm 0.00$
MuVI$0.1$ (3 inf. views) | $\mathbf{0.322\pm 0.00}$ | $\mathbf{0.974\pm 0.01}$ | $\mathbf{0.998\pm 0.00}$ | $\mathbf{0.986\pm 0.01}$

Table 1: Performance comparison on the synthetic data. Average metric scores
across five independent runs $\pm$ standard deviation.  Figure 3: Performance
across multiple prior noise configurations and prior information availability.
Raw results of the unmatched factors (top) versus the results after the
corrected order in a post processing step (bottom). Each boxplot summarizes
five independent runs of the same configuration.

####  4.1.3 Training

The parameter and training settings for the baselines are summarized in the
Appendix A. We perform a grid search on the single hyperparameter of MuVI, the
auxiliary constants $\alpha_{j,k}$ for the factor loadings that are not part
of the prior feature sets. Among $\left\\{0.01,0.03,0.05,0.1\right\\}$,
$\alpha_{j,k}=0.03$ reliably performs the best. We train several models with
different noise configuration and prior information availability. We start
with an uninformed model and iteratively increase the amount of the prior
information by informing a single view, two views and eventually three views.
An important consideration when injecting prior information for a subset of
views is assuming that the features of the rest of the views are not present
in the prior feature sets, and informing them accordingly. In other words, we
choose the same $\alpha_{j,k}$ for all the factor loadings for the uninformed
views that we choose for the factor loadings of the informed views that are
absent in the prior feature sets. The main rationale behind this decision is
encouraging the informed views to learn first, and gradually inform the rest
of the uninformed views implicitly via the structure of the shared factors.
Also, we observed empirically that the uninformed views are unrestricted in
terms of learning and typically converge very early during training and are
unable to “unlearn” a suboptimal structure of the factor loadings.

####  4.1.4 Results

We provide a visual summary of the experiments in Figure 1, and a detailed
overview of the metric scores in Table 1. We first notice a difference in the
quality of the latent representation, measured by the RMSE. MuVI, MOFA and the
autoencoder (AE) perform comparably well. BASS appears to sacrifice
reconstruction quality in favor of additional sparsity, which is supported by
a very high precision and relatively lower recall. Our modelling approach, on
the other hand, strikes a better balance between learning a good latent
representation with low reconstruction error and pruning away the majority of
superfluous connections between the latent factors and the observed features.
In addition, the high recall indicates a better capability of MuVI to retrieve
almost all positive signals, at the cost of some redundant signals. In
addition to the uninformed version of our model, we report the results for
three additional model configurations. The first model simulates a realistic
scenario when the available information is very limited and noisy.
Specifically, we inform only the first view with highly perturbed feature sets
with a noise fraction of $50\%$. Compared to the uninformed MuVI model, we
notice an increase the precision, while the recall remains intact, thereby
further improving the $F_{1}$. In the next two models, we first reduce the
amount of noise to $10\%$, then increase the number of informed views to 3. We
observe additional improvements in all three binary scores, which further
emphasizes the benefit of increasing the quality and the quantity of the prior
information.

Next, we assess the robustness of MuVI against severe noise, and its
capability of utilizing useful prior information to infer a corrected latent
representation. Figure 3 depicts a comprehensive overview of the results
across multiple training settings. In the top row we show the results for the
unmatched factors, whereas in the bottom row the factors underwent a
permutation post training to match the true order. For a fair comparison,
since all baselines are completely unidentifiable and require factor matching,
the scores reported in Table 1 also originate from the bottom row. On the
common horizontal axes, from left to right, we increase the number of views
that are informed starting with a single view in the first column, first and
second view in the middle column, and the first three views in the last
column. Within each subplot, the results are computed and documented
separately for each view. Finally, each group of view-specific boxplots is
ordered by the fraction of prior noise in ascending order. We notice several
attributes of our model. First, MuVI exhibits no difficulties in learning a
good representation from mild to moderate noise in the prior information. As
expected, a severe perturbation reduces the ability of our model to find the
proper ordering of the factors dramatically. However, the bottom row provides
evidence that the model merely loses its ability to be identifiable, since the
quality of the representation remains virtually intact. To further reinforce
the idea of partial identifiability offered by our model, notice the increase
in the $F_{1}$ score when moving along the subplots from left to right, that
is, increasing the number of informed views. The informed factors serve as
anchors in the overall ordering, while at the same time decreasing the number
of possible permutations available for potentially uninformed or view-specific
factors. The effect of the prior information is noticeable not only within
each informed view, but also across the combined observations. For instance,
the uninformed views in the first column, i.e. 1, 2 and 3, fluctuate around an
$F_{1}$ score of 0.65, if we ignore the results with severe noise of $90\%$
and $100\%$. In the second and the third column, the results for the
uninformed views alone improve significantly, which suggests that the shared
factors are capable of correctly transferring their prior information among
the uninformed views.

Finally, we look at the ability of MuVI to communicate the underlying
relationship between the views. The bottom heatmap in Figure 2 corresponds to
the MAP estimates of the factor scales $\delta_{k}^{(m)}$ learned during
inference. The model successfully learns to turn off irrelevant factors and
highlights cross-view relationships by assigning large positive values to
active and shared factors, thereby diminishing the regularization penalty
imposed by the sparse prior.

###  4.2 The Cancer Genome Atlas (TCGA)

Figure 4: Results of the TCGA dataset. A t-SNE embedding space learned from
the latent representation of MuVI (left). Inferred latent factors that are
informed by androgen response and melanogenesis a priori (center). Top 15
features based on the absolute value of their corresponding factor loadings
for the shared androgen response factor among mRNA and RPPA (right).

We investigate a large dataset of comprehensive multi-omic profiling of over
11 thousand samples from 33 cancer types[Tomczak et al., 2015]. Each sample
comprises sets of features of various sizes across four biologically distinct
views: DNA methylation, mRNA expression, microRNA expression and reverse phase
protein array (RPPA). We include a preprocessing step of centering the data
and standardizing each view globally due to the large discrepancy in value
ranges across different views. In addition, samples exhibit partial or
complete missingness in individual views, which we accommodate in our modeling
approach. Prior to training MuVI on the TCGA data, we derive a collection of
gene set annotations (acting as feature sets) from curated public databases
such as MSigDB[Liberzon et al., 2015], Reactome[Fabregat et al., 2018] and
KEGG[Kanehisa and Goto, 2000]. A gene set consists of a group of genes that
are biologically meaningful when co-expressed under certain conditions, or
describe functionally distinct pathways in biological systems. As informative
prior we consider gene set annotations where at least 15 genes were present in
the data, resulting in 360 annotations with a median size of 30. Since gene
sets encompass only genomic and transcriptomic features, we are limited to
informing only the first two views, namely the DNA and the mRNA. We take
advantage of the stochastic variational inference algorithm and propagate
batches of 1,000 samples during training. The algorithm terminates when the
optimization objective no longer improves based on the difference between
consecutive iterations, and after a certain number of patience steps. Training
converges in less than 12 minutes on a single NVIDIA Quadro RTX 5000 GPU with
16 GB memory.

In Figure 4 we display a summary of the results on the TCGA dataset. First, we
assess the quality of the latent representation inferred by our model. We
apply a t-distributed stochastic neighbor embedding (t-SNE) approach to
further compress the latent space into two dimensions for better
visualization[Van der Maaten and Hinton, 2008]. In the left we plot the new
embeddings, colored by their corresponding cancer type. The model is able to
learn a meaningful structure by grouping cancer patients of a similar type
closer together. Next, we explore individual latent factors to gain more
insights into the main drivers of heterogeneity with respect to particular
cancer types. Melanogenesis (K) from the KEGG collection and the androgen
response (H) from the hallmark collection describe general biological
processes. Melanogenesis describes the process of producing the melanin
pigments from melanocytes, commonly found in the epidermis and hair follicles.
Androgen response orchestrates the activation of male hormone receptors, and
plays an important role in the development and progression of prostate
cancer[Fujita and Nonomura, 2019]. In the middle of Figure 4 we illustrate a
scatter plot by mapping our samples onto the two latent axes informed and
described by androgen response and melanogenesis. We highlight three cancer
types with the same coloring as the previous t-SNE plot, while fading the
remaining 30 cancer types to gray. As anticipated, melanogenesis helps
identify the only two groups of samples from tissues with skin cutaneous
melanoma (SKCM) and uveal melanoma (UVM). Likewise, androgen response assigns
significantly higher scores to patients with prostate adenocarcinoma (PRAD).
Next, in the rightmost plots, we dive one level deeper and inspect the top 15
features of the androgen response factor, which is active and shared among the
mRNA and the RPPA view. We decide whether a factor is active in view $m$ if
its corresponding factor scale is not virtually zero, i.e.
$\delta_{k}^{(m)}>0.01$, as well as if the variance explained by the factor
exceeds a threshold of $0.5\%$, based on the coefficient of determination
$R^{2}$. We reiterate that we informed only mRNA among the two views. Hence,
we expect to see annotated genes weighted highly based on the absolute value
of their corresponding factor loading. A powerful feature of MuVI is the
refinement of provided feature sets by adding or removing features in a data-
driven manner. Investigating inferred genes for the Androgen Resonse factor,
we discover several previously known biomarkers of clinical relevance such as
STEAP1, SLC45A3[Perner et al., 2013] and ACPP. The six-transmembrane
epithelial antigen of the prostate 1 (STEAP1) is overexpressed in prostate
cancer and serves as a diagnostic and prognostic biomarker[Khanna et al.,
2021]. Also, it has been shown that levels of ACPP or prostatic specific acid
phosphatase (PSAP), increase with prostate cancer progression[Kong and Byun,
2013]. We now move on to the uninformed view, where we expect the androgen
response factor of mRNA to influence the structure of the factor loadings in
the RPPA, such that the inferred factor reveals additional biomarkers linked
to prostate cancer. Based on the amount of variance explained by this factor
in the RPPA view, and the significance of the shared and informed factor from
the mRNA view, we have strong evidence about the association of the top
proteins with the disease. For example, inferred BCL2 is an oncoprotein that
inhibits the process of apoptosis, or the process of programmed cell death. It
has been demonstrated that BCL2 is a prerequisite for the progression of
cancer cells from an androgen-dependent stage, the earlier stages of the
cancer progression, to an androgen-independent growth stage[Lin et al., 2007].
Moreover, therapeutic approaches that target the checkpoint kinase 1 (CHEK1)
have demonstrated favorable clinical outcomes for prostate cancer
patients[Drápela et al., 2020].

##  5 Conclusion and Discussion

In this contribution, we have addressed the task of modeling multi-view data
with interpretable latent variable models. We use a Bayesian approach with
structured sparsity to encode domain knowledge in the priors of a factor
analysis model. Our model is able to integrate noisy domain expertise and
results in factors that are inherently interpretable via their pre-defined
feature sets. In contrast to other sparse multi-view models, our approach is
able to recover correct signals from noisy feature sets, while maintaining a
low reconstruction error.
We demonstrate in a real-world application that our model is able to infer
clinically meaningful subpopulations of cancer patients.

## References

  * Ainsworth et al. [2018] Samuel K Ainsworth, Nicholas J Foti, Adrian KC Lee, and Emily B Fox.  oi-vae: Output interpretable vaes for nonlinear group factor analysis.  In _International Conference on Machine Learning_ , pages 119–128. PMLR, 2018.
  * Argelaguet et al. [2018] Ricard Argelaguet, Britta Velten, Damien Arnol, Sascha Dietrich, Thorsten Zenz, John C Marioni, Florian Buettner, Wolfgang Huber, and Oliver Stegle.  Multi-omics factor analysis—a framework for unsupervised integration of multi-omics data sets.  _Molecular systems biology_ , 14(6):e8124, 2018\.
  * Armagan et al. [2011] Artin Armagan, Merlise Clyde, and David Dunson.  Generalized beta mixtures of gaussians.  _Advances in neural information processing systems_ , 24, 2011.
  * Barbie et al. [2009] David A Barbie, Pablo Tamayo, Jesse S Boehm, So Young Kim, Susan E Moody, Ian F Dunn, Anna C Schinzel, Peter Sandy, Etienne Meylan, Claudia Scholl, et al.  Systematic rna interference reveals that oncogenic kras-driven cancers require tbk1.  _Nature_ , 462(7269):108–112, 2009.
  * Bernardo et al. [2003] JM Bernardo, MJ Bayarri, JO Berger, AP Dawid, D Heckerman, AFM Smith, and M West.  Bayesian factor regression models in the “large p, small n” paradigm.  _Bayesian statistics_ , 7:733–742, 2003.
  * Bingham et al. [2019] Eli Bingham, Jonathan P Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D Goodman.  Pyro: Deep universal probabilistic programming.  _The Journal of Machine Learning Research_ , 20(1):973–978, 2019.
  * Buettner et al. [2017] Florian Buettner, Naruemon Pratanwanich, Davis J McCarthy, John C Marioni, and Oliver Stegle.  f-sclvm: scalable and versatile factor analysis for single-cell rna-seq.  _Genome biology_ , 18(1):1–13, 2017.
  * Cao and Fleet [2014] Yanshuai Cao and David J Fleet.  Generalized product of experts for automatic and principled fusion of gaussian process predictions.  _arXiv preprint arXiv:1410.7827_ , 2014.
  * Carvalho et al. [2009] Carlos M Carvalho, Nicholas G Polson, and James G Scott.  Handling sparsity via the horseshoe.  In _Artificial Intelligence and Statistics_ , pages 73–80. PMLR, 2009\.
  * Carvalho et al. [2010] Carlos M Carvalho, Nicholas G Polson, and James G Scott.  The horseshoe estimator for sparse signals.  _Biometrika_ , 97(2):465–480, 2010.
  * Damianou et al. [2012] AC Damianou, Carl Henrik Ek, MK Titsias, and ND Lawrence.  Manifold relevance determination.  In _29th International Conference on Machine Learning, ICML 2012, 26 June 2012 through 1 July 2012, Edinburgh_ , pages 145–152, 2012.
  * Drápela et al. [2020] Stanislav Drápela, Prashant Khirsariya, Wytske M van Weerden, Radek Fedr, Tereza Suchánková, Diana Búzová, Jan Červenỳ, Aleš Hampl, Martin Puhr, William R Watson, et al.  The chk1 inhibitor mu380 significantly increases the sensitivity of human docetaxel-resistant prostate cancer cells to gemcitabine through the induction of mitotic catastrophe.  _Molecular oncology_ , 14(10):2487–2503, 2020\.
  * Engelhardt and Stephens [2010] Barbara E Engelhardt and Matthew Stephens.  Analysis of population structure: a unifying framework and novel methods based on sparse factor analysis.  _PLoS genetics_ , 6(9):e1001117, 2010.
  * Fabregat et al. [2018] Antonio Fabregat, Steven Jupe, Lisa Matthews, Konstantinos Sidiropoulos, Marc Gillespie, Phani Garapati, Robin Haw, Bijay Jassal, Florian Korninger, Bruce May, et al.  The reactome pathway knowledgebase.  _Nucleic acids research_ , 46(D1):D649–D655, 2018\.
  * Fujita and Nonomura [2019] Kazutoshi Fujita and Norio Nonomura.  Role of androgen receptor in prostate cancer: a review.  _The world journal of men’s health_ , 37(3):288–295, 2019.
  * Ghosh et al. [2018] Soumya Ghosh, Jiayu Yao, and Finale Doshi-Velez.  Structured variational learning of bayesian neural networks with horseshoe priors.  In _International Conference on Machine Learning_ , pages 1744–1753. PMLR, 2018.
  * Ghosh et al. [2019] Soumya Ghosh, Jiayu Yao, and Finale Doshi-Velez.  Model selection in bayesian neural networks via horseshoe priors.  _J. Mach. Learn. Res._ , 20(182):1–46, 2019.
  * Hoffman et al. [2013] Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley.  Stochastic variational inference.  _Journal of Machine Learning Research_ , 2013.
  * Hotelling [1992] Harold Hotelling.  Relations between two sets of variates.  In _Breakthroughs in statistics_ , pages 162–190. Springer, 1992\.
  * Kanehisa and Goto [2000] Minoru Kanehisa and Susumu Goto.  Kegg: kyoto encyclopedia of genes and genomes.  _Nucleic acids research_ , 28(1):27–30, 2000\.
  * Khanna et al. [2021] Karan Khanna, Nikki Salmond, Kalan S Lynn, Hon S Leong, and Karla C Williams.  Clinical significance of steap1 extracellular vesicles in prostate cancer.  _Prostate cancer and prostatic diseases_ , 24(3):802–811, 2021.
  * Kingma and Welling [2013] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Klami et al. [2013] Arto Klami, Seppo Virtanen, and Samuel Kaski.  Bayesian canonical correlation analysis.  _Journal of Machine Learning Research_ , 14(4), 2013.
  * Klami et al. [2014] Arto Klami, Seppo Virtanen, Eemeli Leppäaho, and Samuel Kaski.  Group factor analysis.  _IEEE transactions on neural networks and learning systems_ , 26(9):2136–2147, 2014.
  * Kong and Byun [2013] Hoon Young Kong and Jonghoe Byun.  Emerging roles of human prostatic acid phosphatase.  _Biomolecules & therapeutics_, 21(1):10, 2013\.
  * Lan et al. [2014] Andrew S Lan, Andrew E Waters, Christoph Studer, and Richard G Baraniuk.  Sparse factor analysis for learning and content analytics.  _The Journal of Machine Learning Research_ , 15(1):1959–2008, 2014.
  * Lee and van der Schaar [2021] Changhee Lee and Mihaela van der Schaar.  A variational information bottleneck approach to multi-omics data integration.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1513–1521. PMLR, 2021.
  * Liberzon et al. [2015] Arthur Liberzon, Chet Birger, Helga Thorvaldsdóttir, Mahmoud Ghandi, Jill P Mesirov, and Pablo Tamayo.  The molecular signatures database hallmark gene set collection.  _Cell systems_ , 1(6):417–425, 2015.
  * Lin et al. [2007] Yuting Lin, Junichi Fukuchi, Richard A Hiipakka, John M Kokontis, and Jialing Xiang.  Up-regulation of bcl-2 is required for the progression of prostate cancer cells from an androgen-dependent to an androgen-independent growth stage.  _Cell research_ , 17(6):531–536, 2007.
  * Louizos et al. [2017] Christos Louizos, Karen Ullrich, and Max Welling.  Bayesian compression for deep learning.  _Advances in neural information processing systems_ , 30, 2017.
  * MacKay et al. [1994] David JC MacKay et al.  Bayesian nonlinear modeling for the prediction competition.  _ASHRAE transactions_ , 100(2):1053–1062, 1994\.
  * Mitchell and Beauchamp [1988] Toby J Mitchell and John J Beauchamp.  Bayesian variable selection in linear regression.  _Journal of the american statistical association_ , 83(404):1023–1032, 1988.
  * Park and Casella [2008] Trevor Park and George Casella.  The bayesian lasso.  _Journal of the American Statistical Association_ , 103(482):681–686, 2008.
  * Perner et al. [2013] Sven Perner, Niels J Rupp, Martin Braun, Mark A Rubin, Holger Moch, Manfred Dietel, Nicolas Wernert, Klaus Jung, Carsten Stephan, and Glen Kristiansen.  Loss of slc45a3 protein (prostein) expression in prostate cancer is associated with slc45a3-erg gene rearrangement and an unfavorable clinical course.  _International journal of cancer_ , 132(4):807–812, 2013.
  * Piironen and Vehtari [2017] Juho Piironen and Aki Vehtari.  Sparsity information and regularization in the horseshoe and other shrinkage priors.  _Electronic Journal of Statistics_ , 11(2):5018–5051, 2017.
  * Polson and Scott [2010] Nicholas G Polson and James G Scott.  Shrink globally, act locally: Sparse bayesian regularization and prediction.  _Bayesian statistics_ , 9(501-538):105, 2010.
  * Ranganath et al. [2014] Rajesh Ranganath, Sean Gerrish, and David Blei.  Black box variational inference.  In _Artificial intelligence and statistics_ , pages 814–822. PMLR, 2014.
  * Ročková and George [2018] Veronika Ročková and Edward I George.  The spike-and-slab lasso.  _Journal of the American Statistical Association_ , 113(521):431–444, 2018.
  * Shapiro [1985] Alexander Shapiro.  Identifiability of factor analysis: Some results and open problems.  _Linear Algebra and its Applications_ , 70:1–7, 1985.
  * Svensson et al. [2020] Valentine Svensson, Adam Gayoso, Nir Yosef, and Lior Pachter.  Interpretable factor models of single-cell rna-seq via variational autoencoders.  _Bioinformatics_ , 36(11):3418–3421, 2020.
  * Thurstone [1931] Louis Leon Thurstone.  Multiple factor analysis.  _Psychological review_ , 38(5):406, 1931.
  * Tibshirani [1996] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 58(1):267–288, 1996.
  * Tomczak et al. [2015] Katarzyna Tomczak, Patrycja Czerwińska, and Maciej Wiznerowicz.  The cancer genome atlas (tcga): an immeasurable source of knowledge.  _Contemporary oncology_ , 19(1A):A68, 2015.
  * Van der Maaten and Hinton [2008] Laurens Van der Maaten and Geoffrey Hinton.  Visualizing data using t-sne.  _Journal of machine learning research_ , 9(11), 2008.
  * Zhao et al. [2016] Shiwen Zhao, Chuan Gao, Sayan Mukherjee, and Barbara E Engelhardt.  Bayesian group factor analysis with structured sparsity.  _The Journal of Machine Learning Research_ , 2016.

##  Appendix A Appendix A

The parameter and training settings for the baselines are provided here. We
train MOFA combining the ARD and the spike-and-slab prior on the weights for
column-wise and element-wise sparsity and apply BASS with default parameters
and $20$ parameter-expanded expectation maximization (PX-EM) iterations. The
autoencoder consists of an encoder-decoder pair for each observed view,
connected via a common latent space governed by a product of experts (PoE)
approach[Lee and van der Schaar, 2021, Cao and Fleet, 2014]. The encoders
comprise two hidden layers of size 64 and 32, each followed by a ReLU
activation, whereas the decoders perform a linear mapping of the latent code
to the observed features of each view. To achieve sparsity in the factor-to-
feature mapping, we introduce an L1 penalty term in the respective decoder
weights.

##  Appendix B Appendix B

We assessed the effect of the threshold on precision and recall by generating
PR curves. Figure 5 depicts the precision-recall curves for all models.

Figure 5: Precision-recall curves for MuVI0.5 and the baseline models: BASS,
MOFA, AE.

##  Appendix C Appendix C

Here, we show that our approach is robust to the choice of the single
hyperparameter $\alpha$. Figure 6 compiles the results of a hyperparameter
sensitive analysis for the synthetic experiments in Section 4.1.

Figure 6: Sensitivity analysis to the auxiliary $\alpha_{j,k}$ set to inactive
signals in the prior feature set. Results when informing only the first view.

##  Appendix D Appendix D

Finally, we show in Figure 7 that an autoencoder does not capture any
additional information than MuVI in the latent space in terms of identifying
subpopulations.

Figure 7: Embedding space learned by the hybrid autoencoder.