# Multi-study factor regression model: an application in nutritional
epidemiology

Roberta De Vito
Dept. of Biostatistics  Acknowledgments: RDV was supported by the US National
Institutes of Health, NIGMS/NIH COBRE CBHD P20GM109035. This Manuscript was
prepared using HCHSSOL Research Materials obtained from the NHLBI Biologic
Specimen and Data Repository Information Coordinating Center and does not
necessarily reflect the opinions or views of the HCHSSOL or the NHLBI. The
authors would like to thank David Rossell and Giovanni Parmigiani for their
precious advice.  School of Public Health  Brown University  Providence  RI
Data Science Initiative  Brown University  Providence  RI
and
Alejandra Avalos-Pacheco
Applied Statistics Research Unit  TU Wien  Vienna  Austria
Harvard-MIT Center for Regulatory Science  Harvard University  Boston  MA

###### Abstract

Diet is a risk factor for many diseases. In nutritional epidemiology, studying
reproducible dietary patterns is critical to reveal important associations
with health. However, it is challenging: diverse cultural and ethnic
backgrounds may critically impact eating patterns, showing heterogeneity,
leading to incorrect dietary patterns and obscuring the components shared
across different groups or populations.

Moreover, covariate effects generated from observed variables, such as
demographics and other confounders, can further bias these dietary patterns.
Identifying the shared and group-specific dietary components and covariate
effects is essential to drive accurate conclusions.

To address these issues, we introduce a new modeling factor regression, the
Multi-Study Factor Regression (MSFR) model. The MSFR model analyzes different
populations simultaneously, achieving three goals: capturing shared
component(s) across populations, identifying group-specific structures, and
correcting for covariate effects.

We use this novel method to derive common and ethnic-specific dietary patterns
in a multi-center epidemiological study in Hispanic/Latinos community. Our
model improves the accuracy of common and group dietary signals and yields
better prediction than other techniques, revealing significant associations
with health.

In summary, we provide a tool to integrate different groups, giving accurate
dietary signals crucial to inform public health policy.

Keywords: Factor analysis, Joint analysis, ECM algorithm, Nutritional
epidemiology, Dietary patterns

##  1 Introduction

The impact of diet on diseases such as cancer (Willett and MacMahon, 1984),
obesity (Foster et al., 2003) and cardiovascular disease (Kromhout, 2001; GBD
2013 Risk Factors Collaborators and others, 2015), has continually affected
the world, including the United States, accounting for over 500,000 deaths
annually, with 280,000 of those deaths due to obesity (Aune et al., 2016).
Therefore, detecting eating patterns is crucial to assess the association
between diet and diseases to guide and inform public health policy.

Diverse cultural and ethnic background components may have crucial impacts on
eating patterns. In nutritional epidemiology, heterogeneous high-dimensional
data coming from different populations or ethnic backgrounds (Stearns et al.,
2017; Zulyniak et al., 2017) often arise (Haines et al., 1999), providing an
opportunity to estimate shared dietary patterns and differences across them,
leading to a better understanding of their association with certain types of
diseases.

When analyzing different populations or groups, some diets and foods are
reproducible and shared, while other diet characteristics are specific to each
population. To handle this task, standard procedures involve stacking all the
populations and performing traditional techniques such as factor analysis
(Haines et al., 1999; Schulze et al., 2003; Edefonti et al., 2012; Bennett et
al., 2022) to estimate dietary patterns, or performing factor analysis in each
population and then adopting a measure of correlation to identify the most
similar dietary patterns across the different populations or subgroups
(Castello et al., 2014; Castelló et al., 2022). However, identifying dietary
patterns in different populations can be challenging, mainly due to the
diversity of ethnicities and demographic compositions among populations.
Consumption patterns of foods, recipes or beverages can vary significantly
between populations, critically affecting the risk of disease outcomes.
Population-based identification of these diets is crucial to provide essential
knowledge and information for broad public health policies. However,
distinctly unique populations can impact and mask the shared and reproducible
dietary patterns among all the populations. Therefore, there is a crucial need
for rigorous statistical techniques to identify concurrently common and study-
specific diet characteristics to clearly estimate the association between diet
and diseases.

Currently, researchers in the field are focused on identifying group or
reproducible signal patterns. A recent development in this direction is the
multi-study factor analysis (MSFA) model (De Vito et al., 2019), which is able
to identify common and study-specific characteristics in gene expression data
across multiple studies.

This approach was also adopted in nutritional epidemiology. For instance, De
Vito et al. (2019) analyzed the International Head and Neck Cancer
Epidemiology (INHANCE) consortium (Hashibe et al., 2007; Conway et al., 2009)
data via MSFA, revealing important common and population-specific diet
characteristics and their association with the two cancer types. In this
direction, De Vito et al. (2022) adopted a multi-study Bayesian approach in a
nutritional epidemiological setting in 4 US field sites (Bronx, Chicago,
Miami, San Diego) in the Hispanic Community Health Study/Study of Latinos
(HCHS/SOL) (Sorlie et al., 2010). More recently, Stephenson et al. (2020)
introduced a generalization of mixture models, namely the Robust Profile
Clustering (RPC), to handle diverse populations or studies in a case-control
birth defects study in the United States (Yoon et al., 2001). The RPC
identifies robust food “global” clusters across all populations and locally
population clusters.

All of these approaches highlight the crucial need in nutritional epidemiology
to identify common and study diet components concurrently, in order to gain a
better understanding of their association with disease outcomes. However, none
of these approaches accounts for covariates or confounders essential to
identify the diet signals accurately. Specifically, the diet is affected by
numerous confounders (Trichopoulos et al., 1985; Tucker et al., 1997), such as
smoking, alcohol, gender, and other covariates, which can obscure both the
common and study-specific signals. Thus, not adjusting for multiple potential
confounding factors could result in inaccurate conclusions.

In the standard factor analysis framework, models that account for observed
covariates are well-studied. For instance, West, 2003 introduced sparse latent
factor models with many explanatory variables. This approach enables a more
precise assessment of the factor structure, leading to improved breast cancer
data predictions. Carvalho et al., 2008 developed a factor regression model
that estimated both factor loadings and regression coefficients. This paper
provided a modeling approach that improved the characterization of the factor
loading patterns in cancer genomic applications by accounting for observed
covariates. This approach also showed an increased predictive ability when
considering observed covariates. More recently, Avalos-Pacheco et al., 2022
introduced the Bayesian factor regression model, a novel approach that jointly
analyzes multiple studies and adjusts for systematic biases, such as batch
effects, in a latent factor regression framework. The authors illustrated
their analysis’ power in both signal detection and predictive ability and
empirically showed that inaccurate conclusions can be drawn if batch effects
are not appropriately considered.

Although all these methods revealed the crucial need to account for both
latent structures and observed variables, none of these methods considered
study-specific latent patterns that contribute significantly to the variations
of each population. Indeed, these populations or background dietary patterns
can reveal essential differences among groups and sometimes mask the shared
signal among groups, drastically biasing the associations with the disease.
Therefore, it is critical to identify both the shared and study-specific
latent components while accounting for other observed covariates, like
counfonders and/or batch effects. To address this challenge, we propose a new
modeling approach, the Multi-Study Factor Regression (MSFR) model, which is
able to analyze multiple studies/ populations/groups jointly and estimate
three different components: (1) the reproducible and shared latent components,
(2) the study-specific latent component for each population, (3) the
regression coefficients, which are key to detect confounders.

Our work focuses on the Hispanic Community Health Study/Study of Latinos
(HCHS/SOL) (National Heart, Lung, and Blood Institute and others, 2009), which
is a multi-center epidemiologic study in Hispanic/Latino populations. The
HCHS/SOL aims to detect crucial factors impacting the health of
Hispanics/Latinos (Sorlie et al., 2010). One of the study’s primary goal is to
determine how diet can impact cardiovascular disease (Daviglus et al., 2014)
in different Hispanic/Latino origins, specifically in the Cuban, Puerto Rican,
Dominican, Mexican, and Central/South American backgrounds. Our MSFR model
aims to identify the reproducible dietary components among the different
Hispanic/Latino backgrounds and background dietary components while accounting
for confounders (i.e., smoking, alcohol, etc.) that can obscure the two
dietary signals. In this way, we can assess a more accurate association
between dietary patterns and cardiovascular diseases, such as obesity.

The plan of the paper is as follows. Section 2 motivates our new approach by
describing the Hispanic Community Health Study/Study of Latinos (HCHS/SOL) and
all the variables considered in our analysis. Section 3 presents the multi-
study regression approach, including the model definition and model
estimation. Section 4 presents the performance of our MSFR in different
simulation settings, comparing it with the multi-study factor model, the
factor regression model and a two-step approach aimed to correct for covariate
effects. Section 5 includes the application of MSFR in the HCHS/SOL study,
which motivated our novel method. Section 6 includes the discussion.

##  2 The Hispanic Community Health Study / Study of Latinos

###  2.1 The study

The Hispanic Community Health Study/Study of Latinos (HCHS/SOL) is the largest
multi-site cohort study of Hispanics/Latinos in the United States. The aim of
the study is to provide previously unavailable information on the health
status and disease burden of US Hispanics/Latinos and to explore the
associations between risk factors and disease incidence during the follow-up.
In this manuscript, we consider food and diet as a potential risk for health
index, particularly the healthy eating index and cardiovascular disease. The
study collected 16,415 adults aged 18-74 years residing in four different
sites (Bronx, Chicago, Miami, San Diego) from seven Hispanic/Latino ethnic
backgrounds (Cuban, Dominican Republic, Mexican, Puerto Rican, Central, South
American, and others). The information is collected at baseline (2008-2011)
using a stratified 2-stage probability sampling design, described in detail by
LaVange et al. (2010).

The cohort study collects data via two 24-hour recalls, conducted in person at
the baseline and via telephone, $\leq 30$ days from baseline in the
participant’s preferred language (English 20% and Spanish 80% of the
participants) by trained interviewers. The data is obtained after the two
recalls, and based on thousands of common and Hispanic/Latino foods, contains
a total of 139 nutrients, nutrient ratios, food-group serving counts, and
other food components. All participants (99%) provided at least one recall.

###  2.2 Selection of subjects and data preprocessing

A total of 14,002 participants’ data were collected. We excluded participants
who self-identified as belonging to other ethnic backgrounds and/or present
recall data unreliable (i.e., vitamin E $<0$) and/or show inconsistent (i.e.,
$<1st\mbox{ or }>99th$percentiles) estimates of energy intake. We also removed
individuals with missing ethnic background information. After applying all of
the above-mentioned exclusion criteria, we performed our analyses on a total
of 12,977 participants.

Following De Vito et al. (2022), we selected 53 nutrients that best represent
the overall diet for Hispanics/Latinos. We expanded the fat profile for better
identifying cardiovascular-disease-related dietary habits (Islam et al., 2019;
Lichtenstein, 2003). In particular, we constructed four variables: 1- SCSFA
short-chain saturated fatty acids (i.e., butyric acid), 2- MCSFA medium-chain
saturated fatty acids (i.e., the sum of caproic acid, caprylic acid, capric
acid, and lauric acid), 3- LCSFA long-chain saturated fatty acids (i.e., the
sum of myristic acid, palmitic acid, margaric acid, stearic acid and arachidic
acid), 4- LCSFA long-chain monounsaturated fatty acids (i.e., the sum of
myristoleic acid, palmitoleic acid, oleic acid, and gadoleic acid).

We then derived nutrient intakes from the mean of two available reliable
recalls for including all the diet information. We log-transformed the
nutrient intakes and standardized the data.

An overview of socio-demographic characteristics for each hispanic background
group is provided in Table 1. Note that we only report socio-demographic
characteristics included in our model (Section 5).

| Dominican |  Central American | Cuban | Mexican |  Puerto Rican |  South American
---|---|---|---|---|---|---
| ($n=1257$) | ($n=1443$) | ($n=2126$) | ($n=4940$) | ($n=2314$) | ($n=897$)
Gender |  |  |  |  |  |
Female | 817 | 862 | 1139 | 3031 | 1350 | 530
Male | 440 | 583 | 987 | 1909 | 964 | 367
Education |  |  |  |  |  |
$<$ High School | 538 | 595 | 486 | 2147 | 903 | 212
= High School | 249 | 332 | 627 | 1300 | 609 | 232
$>$ High School | 470 | 518 | 1013 | 1493 | 802 | 453
Alcohol Use |  |  |  |  |  |
Never | 619 | 858 | 1159 | 2521 | 1226 | 461
Former | 575 | 521 | 843 | 2164 | 952 | 410
Current | 63 | 66 | 124 | 255 | 136 | 26
Cigarette |  |  |  |  |  |
0 years | 947 | 999 | 1077 | 3241 | 1122 | 588
$>0$ and $<10$ years | 170 | 294 | 334 | 1240 | 541 | 196
10+ years | 140 | 152 | 715 | 459 | 651 | 113
BMI∗ mean(sd) | 29.4(5.7) | 29.7(5.6) | 29.2(5.7) | 29.8(5.7) | 30.9(6.7) | 28.7 (5.2)
Table 1: Social and Demographic characteristic of the HCHS/SOL consortium for
each hispanic background. ∗ Body Mass Index (kg/$m^{2}$)

Hispanic/Latino backgrounds were not equally distributed at each study site.
Specifically, the population with Mexican (South American) background is the
largest (smaller) among the groups. The majority of individuals were female,
never alcohol use, smoking no cigarettes. The category of Education most
concentrated was “less than high school”, except for the Cuban background with
a concentration on “greater than high school”.

We consider a total of 53 food variables and 5 confounders to best identify
the common signals across the 6 different ethnic background groups and the
study-specific components for each of these groups.

##  3 Multi-study factor regression model

We jointly model the data from $S$ different groups, populations, or studies
that observe the same $p$ responses. We consider vectors
$\mathbf{x}_{is}\in\mathbb{R}^{p}$, observed for $i=1,\dots,n$ individuals in
$s=1,\dots,S$ studies. Our method generalizes the traditional multi-study
factor analysis (MSFA) (De Vito et al., 2021).

Traditional MSFA models $\mathbf{x}_{is}$ as a linear combination of low-
dimensional common factors $\mathbf{f}_{is}\in\mathbb{R}^{q}$, and group-
specific factors $\mathbf{l}_{is}\in\mathbb{R}^{q_{s}}$, $q+q_{s}\ll p$:

$\displaystyle\mathbf{x}_{is}=\boldsymbol{\Phi}\mathbf{f}_{is}+\boldsymbol{\Lambda}_{s}\mathbf{l}_{is}+\mathbf{e}_{is},$
(1)

where $\boldsymbol{\Lambda}_{s}\in\mathbb{R}^{p\times p_{s}}$ for
$s=1,\dots,S$, is the specific loadings matrix,
$\boldsymbol{\Phi}\in\mathbb{R}^{p\times q}$ is the factor loadings matrix,
and $\mathbf{e}_{is}\in\mathbb{R}^{p}$ is the error, distributed as
$\mathbf{e}_{is}\sim N(0,\boldsymbol{\Psi}_{s}),s=1,\dots,S$ where
$\boldsymbol{\Psi}_{s}=\text{diag}\\{\psi_{js}\\}_{j=1}^{p}$ is the
idiosyncratic precision matrix of study $s$. Common and study-specific factors
are assumed to be standard normal and independent between them and the errors.
Marginally, observations $\mathbf{x}_{is}$ are normally distributed, centered
at zero and with covariance matrix
$\boldsymbol{\Sigma}_{s}=\boldsymbol{\Phi}\boldsymbol{\Phi}^{\top}+\boldsymbol{\Lambda}_{s}\boldsymbol{\Lambda}_{s}^{\top}+\boldsymbol{\Psi}_{s}$.
The covariance matrix allows decomposing the total variance into three
distinct parts, i.e. the variance of the common factors
$\boldsymbol{\Phi}\boldsymbol{\Phi}^{\top}$, of the study-specific factors
$\boldsymbol{\Lambda}_{s}\boldsymbol{\Lambda}_{s}^{\top}$ and the error
$\boldsymbol{\Psi}_{s}$.

Model (1) assumes that the data is explained by unobserved latent factors
only, without taking into account prespecified variables that can affect the
outcome in practice. For instance, in our HCHS/SOL application, the MSFA
ignores the effect that observed covariates, such as BMI, alcohol and
cigarette use, can have on cardiovascular-disease-related dietary habits. In
other words, MSFA models the data without allowing for other additional
regression effects that do not affect the underlying dietary patterns given in
the common and study-specific factors. We will later describe how much
inference can suffer when not incorporating such covariate adjustments (Tables
2 and 3, Figures S1, S2, S3, and S4, Supplementary Materials). We tackle this
issue by extending Model (1) to incorporate covariate adjustments by adopting
a factor regression approach (Avalos-Pacheco et al., 2022; Avalos-Pacheco,
2018; Carvalho et al., 2008). Let $\mathbf{b}_{is}$ be $p_{b}$-dimensional
vector of known covariates. The Multy-study Factor Regression model is:

$\mathbf{x}_{is}=\boldsymbol{\beta}\mathbf{b}_{is}+\boldsymbol{\Phi}\mathbf{f}_{is}+\boldsymbol{\Lambda}_{s}\mathbf{l}_{is}+\mathbf{e}_{is}$
(2)

where $\boldsymbol{\beta}\in\mathbb{R}^{p\times p_{b}}$ is the matrix of
regression coefficients.

We remark that the utility of this model goes beyond our Hispanic Community
Health Study. For instance, this model can easily adjust for several covariate
effects, such as complex batch effects, by simply setting
$\mathbf{b}_{is}\in[0,1]^{S}$ as a study indicator. In this setting, $\beta$
captures the mean batch effect adjustment and $\boldsymbol{\Psi}_{s}$ the
variance adjustment.

## Maximum likelihood estimation by ECM algorithm

We fit Model (2) using a deterministic Expectation/Conditional Maximization
(ECM) optimization (Meng and Rubin, 1993), along the lines of De Vito et al.
(2019) to estimate the parameters
$\boldsymbol{\theta}=\left\\{\boldsymbol{\Phi},\boldsymbol{\Lambda}_{s},\boldsymbol{\Psi}_{s},{\color[rgb]{0,0,0}\boldsymbol{\beta}}\right\\}$.
Our proposed ECM, summarized in Algorithm 1, maximises the expected complete-
data log-likelihood
$l_{c}(\boldsymbol{\theta})=E_{\mathbf{F},\mathbf{L}_{s}\mid\boldsymbol{\theta},\mathbf{X}}[p(\boldsymbol{\theta}\mid\mathbf{F},\mathbf{L}_{s},\mathbf{X})]$
iteratively. For simplicity, let
${\widetilde{\mathbf{x}}_{is}}=\mathbf{x}_{is}{\color[rgb]{0,0,0}-\boldsymbol{\beta}\mathbf{b}_{is}}$,
and
$\widehat{\boldsymbol{\theta}}=\left\\{\widehat{\boldsymbol{\Phi}},\widehat{\boldsymbol{\Lambda}}_{s},\widehat{\boldsymbol{\Psi}}_{s},\widehat{\boldsymbol{\beta}}\right\\}$
be the current value of the estimated parameters.

initialize
$\widehat{\boldsymbol{\Phi}}=\boldsymbol{\Phi}^{(0)},\widehat{\boldsymbol{\Lambda}}_{s}=\boldsymbol{\Lambda}_{s}^{(0)},\widehat{\boldsymbol{\Psi}}_{s}=\boldsymbol{\Psi}_{s}^{(0)},\widehat{\boldsymbol{\beta}}=\boldsymbol{\beta}^{(0)}$

while  _$t <T$ and $\epsilon>\epsilon^{*}$_ do

E-step§:

$\begin{array}[]{cccc}&E[\mathbf{f}_{is}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}]=\mathbf{E}_{f_{i}}=\boldsymbol{\delta}{\widetilde{\mathbf{x}}_{is}},&&E[\mathbf{l}_{is}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}]=\mathbf{E}_{l_{is}}=\boldsymbol{\delta}_{s}{\widetilde{\mathbf{x}}_{is}},\\\
&\rm{E}\left[\mathbf{C}_{\widetilde{x}_{s}\widetilde{x}_{s}}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}\right]=\mathbf{E}_{\widetilde{x}_{s}\widetilde{x}_{s}}=\mathbf{C}_{\widetilde{x}_{s}\widetilde{x}_{s}},&&\rm{E}\left[\mathbf{C}_{\widetilde{x}_{s}l_{s}}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}\right]=\mathbf{E}_{\widetilde{x}_{s}l_{s}}=\mathbf{C}_{\widetilde{x}_{s}\widetilde{x}_{s}}\boldsymbol{\delta}_{s}^{\top},\\\
&\rm{E}\left[\mathbf{C}_{\widetilde{x}_{s}f}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}\right]=\mathbf{E}_{\widetilde{x}_{s}f}=\mathbf{C}_{\widetilde{x}_{s}\widetilde{x}_{s}}\boldsymbol{\delta}^{\top},&&\rm{E}\left[\mathbf{C}_{ff}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}\right]=\mathbf{E}_{ff}=\boldsymbol{\delta}\mathbf{C}_{\widetilde{x}_{s}\widetilde{x}_{s}}\boldsymbol{\delta}^{\top}+\boldsymbol{\Delta},\\\
&\rm{E}\left[\mathbf{C}_{l_{s}l_{s}}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}\right]=\mathbf{E}_{l_{s}l_{s}}=\boldsymbol{\delta}_{s}\mathbf{C}_{\widetilde{x}_{s}\widetilde{x}_{s}}\boldsymbol{\delta}_{s}^{\top}+\boldsymbol{\Delta}_{s}&&\rm{E}\left[\mathbf{C}_{fl_{s}}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}\right]=\mathbf{E}_{fl_{s}}=\boldsymbol{\delta}\mathbf{C}_{\widetilde{x}_{s}\widetilde{x}_{s}}\boldsymbol{\delta}_{s}^{\top}+\boldsymbol{\Delta}_{fl}.\end{array}$

CM-step:

$\widehat{\boldsymbol{\Psi}}_{s}=$ |  diag$\left\\{\mathbf{E}_{\widetilde{x}_{s}\widetilde{x}_{s}}+\widehat{\boldsymbol{\Phi}}\mathbf{E}_{ff}\widehat{\boldsymbol{\Phi}}^{\top}+\widehat{\boldsymbol{\Lambda}}_{s}\mathbf{E}_{l_{s}l_{s}}\widehat{\boldsymbol{\Lambda}}_{s}^{\top}-2\mathbf{E}_{\widetilde{x}_{s}f}\widehat{\boldsymbol{\Phi}}^{\top}-\Biggl{.}2\mathbf{E}_{\widetilde{x}_{s}l_{s}}\widehat{\boldsymbol{\Lambda}}_{s}^{\top}+2\widehat{\boldsymbol{\Phi}}\mathbf{E}_{fl_{s}}\widehat{\boldsymbol{\Lambda}}_{s}^{\top}\right\\}$
---|---
$\mbox{vec}(\widehat{\boldsymbol{\Phi}})=$ | $\sum_{s=1}^{S}\left(\mathbf{E}_{ff}^{\top}\otimes n_{s}\widehat{\boldsymbol{\Psi}}^{-1}_{s}\right)\mbox{vec}\left(n_{s}\widehat{\boldsymbol{\Psi}}^{-1}_{s}\mathbf{E}_{\widetilde{x}_{s}f}-n_{s}\widehat{\boldsymbol{\Psi}}^{-1}_{s}\widehat{\boldsymbol{\Lambda}}_{s}\mathbf{E}_{fl_{s}}^{\top}\right)$
$\widehat{\boldsymbol{\Lambda}}_{s}=$ | $\Biggl{(}\mathbf{E}_{x_{s}f}-\widehat{\boldsymbol{\Phi}}\mathbf{E}_{fl_{s}}\Biggr{)}\Biggl{(}\mathbf{E}_{l_{s}l_{s}}\Biggr{)}^{-1}$
$\widehat{\boldsymbol{\beta}}=$ | $\left[\sum_{s=1}^{S}\sum_{i}^{n_{s}}(\mathbf{x}_{is}-\widehat{\boldsymbol{\Phi}}\mathbf{E}_{f_{i}}-\widehat{\boldsymbol{\Lambda}}_{s}\mathbf{E}_{l_{is})}\mathbf{b}_{is}^{\top}\right]\left[\sum_{s=1}^{S}\sum_{i}^{n_{s}}(\mathbf{b}_{is}\mathbf{b}_{is}^{\top})\right]^{-1}$

set $\boldsymbol{\theta}^{(t+1)}=\widehat{\boldsymbol{\theta}}$

compute $t=t+1$ and $\epsilon=\mid l^{(t+1)}-l^{(t)}\mid$, with
$l^{(t+1)}=l_{c}(\boldsymbol{\theta}^{(t)})+\frac{l_{c}(\boldsymbol{\theta}^{(t+1)})-l_{c}(\boldsymbol{\theta}^{(t)})}{1-c^{(t)}}$,
and
$c^{(t)}=\frac{l_{c}(\boldsymbol{\theta}^{(t+1)})-l_{c}(\boldsymbol{\theta}^{(t)})}{l_{c}(\boldsymbol{\theta}^{(t)})-l_{c}(\boldsymbol{\theta}^{(t-1)})}$

end while

§ where
$\begin{array}[]{cccccc}&\mathbf{C}_{\widetilde{x}_{s}\widetilde{x}_{s}}=\dfrac{\sum_{i=1}^{n_{s}}{\widetilde{\mathbf{x}}_{is}}{\widetilde{\mathbf{x}}_{is}}^{\top}}{n_{s}},&&\mathbf{C}_{ff}=\dfrac{\sum_{i=1}^{n_{s}}\mathbf{f}_{is}\mathbf{f}_{is}^{\top}}{n_{s}},&&\mathbf{C}_{l_{s}l_{s}}=\dfrac{\sum_{i=1}^{n_{s}}\mathbf{l}_{is}\mathbf{l}_{is}^{\top}}{n_{s}},\\\
&\mathbf{C}_{\widetilde{x}_{s}f}=\dfrac{\sum_{i=1}^{n_{s}}{\widetilde{\mathbf{x}}_{is}}\mathbf{f}_{is}^{\top}}{n_{s}},&&\mathbf{C}_{\widetilde{x}_{s}l_{s}}=\dfrac{\sum_{i=1}^{n_{s}}{\widetilde{\mathbf{x}}_{is}}\mathbf{l}_{is}^{\top}}{n_{s}},&&\mathbf{C}_{fl_{s}}=\dfrac{\sum_{i=1}^{n_{s}}\mathbf{f}_{is}\mathbf{l}_{is}^{\top}}{n_{s}},\end{array}$

$\boldsymbol{\Delta}_{fl}=\rm{Cov}[\mathbf{l}_{is},\mathbf{f}_{is}|{\widetilde{\mathbf{x}}_{is}},\boldsymbol{\theta}_{t-1}]=\boldsymbol{\Phi}^{\top}\boldsymbol{\Sigma}_{s}^{-1}\boldsymbol{\Lambda}_{s},$
and
$\begin{array}[]{cccc}&\boldsymbol{\delta}=\boldsymbol{\Phi}^{\top}\boldsymbol{\Sigma}_{s}^{-1},&&\boldsymbol{\Delta}=\rm{Var}[\mathbf{f}_{is}|{\widetilde{\mathbf{x}}_{is}},\boldsymbol{\theta}_{t-1}]=\mathbf{I}_{q}-\boldsymbol{\Phi}^{\top}\boldsymbol{\Sigma}_{s}^{-1}\boldsymbol{\Phi},\\\
&\boldsymbol{\delta}_{s}=\boldsymbol{\Lambda}_{s}^{\top}\boldsymbol{\Sigma}_{s}^{-1},&&\boldsymbol{\Delta}_{s}=\rm{Var}[\mathbf{l}_{is}|{\widetilde{\mathbf{x}}_{is}},\boldsymbol{\theta}_{t-1}]=\mathbf{I}_{q_{s}}-\boldsymbol{\Lambda}_{s}^{\top}\boldsymbol{\Sigma}_{s}^{-1}\boldsymbol{\Lambda}_{s}.\end{array}$

Algorithm 1 ECM algorithm for Multi-study Factor Regression

In the E-step, we obtained the expected complete-data log-likelihood

$\begin{split}l_{c}(\boldsymbol{\theta})=&\sum_{s=1}^{S}\Biggl{\\{}-\frac{n_{s}}{2}\log|\boldsymbol{\Psi}_{s}|-\frac{n_{s}}{2}tr\left[\boldsymbol{\Psi}_{s}^{-1}\frac{\sum_{i}^{n_{s}}\left(E[{\widetilde{\mathbf{x}}_{is}}{\widetilde{\mathbf{x}}_{is}}^{\top}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}]\right.}{n_{s}}\right.\Biggr{.}\\\
&\frac{+\boldsymbol{\Phi}E[\mathbf{f}_{is}\mathbf{f}_{is}^{\top}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}]\boldsymbol{\Phi}^{\top}+\boldsymbol{\Lambda}_{s}E[\mathbf{l}_{is}\mathbf{l}_{is}^{\top}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}]\boldsymbol{\Lambda}_{s}^{\top}-2E[{\widetilde{\mathbf{x}}_{is}}\mathbf{f}_{is}^{\top}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}]\boldsymbol{\Phi}^{\top})}{n_{s}}\\\
&\Biggl{.}\frac{-2E[{\widetilde{\mathbf{x}}_{is}}\mathbf{l}_{is}^{\top}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}]\boldsymbol{\Lambda}_{s}^{\top}+2\boldsymbol{\Phi}E[\mathbf{f}_{is}\mathbf{l}_{is}^{\top}|{\widetilde{\mathbf{x}}_{is}},\widehat{\boldsymbol{\theta}}]\boldsymbol{\Lambda}_{s}^{\top})}{n_{s}}\Biggr{]}\Biggr{\\}}+C,\end{split}$
(3)

where $C$ is a constant. Equation (3) depends on the latent factors
$\mathbf{f}_{is}$ and $\mathbf{l}_{is}$ in terms of their conditional mean and
their conditional second moments, shown in Algorithm 1 and Supplement A. All
the moments are computed based on the multivariate normal distribution.

Note that the conditional expectations require to invert the $p\times p$
covariate matrix $\widehat{\boldsymbol{\Sigma}}_{s}$, through

$\begin{array}[]{cccc}&\boldsymbol{\delta}=\widehat{\boldsymbol{\Phi}}^{\top}\widehat{\boldsymbol{\Sigma}}_{s}^{-1},&&\boldsymbol{\delta}_{s}=\widehat{\boldsymbol{\Lambda}}_{s}^{\top}\widehat{\boldsymbol{\Sigma}}_{s}^{-1},\end{array}$
(4)

which can be cumbersome and challenging as the number of the variables $p$
increases. To address this issue, we apply the Woodbury matrix identity
(Morgenstern and Woodbury, 1950) twice, obtaining

$\begin{array}[]{cccc}&\boldsymbol{\delta}=(\mathbf{I}_{q}+\widehat{\boldsymbol{\Phi}}^{\top}\widehat{\mathbf{W}}_{\boldsymbol{\Lambda}_{s}}\widehat{\boldsymbol{\Phi}})^{-1}\widehat{\boldsymbol{\Phi}}^{\top}\widehat{\mathbf{W}}_{\boldsymbol{\Lambda}_{s}},&&\boldsymbol{\delta}_{s}=(\mathbf{I}_{q_{s}}+\widehat{\boldsymbol{\Lambda}}_{s}^{\top}\widehat{\mathbf{W}}_{\boldsymbol{\Phi}}\widehat{\boldsymbol{\Lambda}}_{s})^{-1}\widehat{\boldsymbol{\Lambda}}_{s}^{\top}\widehat{\mathbf{W}}_{\boldsymbol{\Phi}},\end{array}$
(5)

with

$\begin{split}\widehat{\mathbf{W}}_{\boldsymbol{\Lambda}_{s}}=(\widehat{\boldsymbol{\Lambda}}_{s}\widehat{\boldsymbol{\Lambda}}_{s}^{\top}+\widehat{\boldsymbol{\Psi}}_{s})^{-1}=&\widehat{\boldsymbol{\Psi}}_{s}^{-1}-\widehat{\boldsymbol{\Psi}}_{s}^{-1}\widehat{\boldsymbol{\Lambda}}_{s}(\mathbf{I}_{q_{s}}+\widehat{\boldsymbol{\Lambda}}_{s}^{\top}\widehat{\boldsymbol{\Psi}}_{s}^{-1}\widehat{\boldsymbol{\Lambda}}_{s})^{-1}\widehat{\boldsymbol{\Lambda}}_{s}^{\top}\widehat{\boldsymbol{\Psi}}_{s}^{-1}\\\
\widehat{\mathbf{W}}_{\boldsymbol{\Phi}}=(\widehat{\boldsymbol{\Phi}}\widehat{\boldsymbol{\Phi}}^{\top}+\widehat{\boldsymbol{\Psi}}_{s})^{-1}=&\widehat{\boldsymbol{\Psi}}_{s}^{-1}-\widehat{\boldsymbol{\Psi}}_{s}^{-1}\widehat{\boldsymbol{\Phi}}(\mathbf{I}_{q}+\widehat{\boldsymbol{\Phi}}^{\top}\widehat{\boldsymbol{\Psi}}_{s}^{-1}\widehat{\boldsymbol{\Phi}})^{-1}\widehat{\boldsymbol{\Phi}}^{\top}\widehat{\boldsymbol{\Psi}}_{s}^{-1}.\end{split}$
(6)

We have carried the inference of $\mathbf{E}_{f_{is}}$ and
$\mathbf{E}_{l_{is}}$ using Equation (6).

In the CM-step, we maximize $l_{c}(\boldsymbol{\theta})$ with respect to
$\boldsymbol{\Phi},\boldsymbol{\Lambda}_{s},\boldsymbol{\Psi}_{s},{\color[rgb]{0,0,0}\boldsymbol{\beta}}$.
For simplicity, we can re-write Equation (3) as

$\begin{split}l_{c}(\boldsymbol{\theta})=&\sum_{s=1}^{S}\left(-\frac{n_{s}}{2}\log|\boldsymbol{\Psi}_{s}|-\frac{n_{s}}{2}tr\left[\boldsymbol{\Psi}_{s}^{-1}\left(\mathbf{E}_{\widetilde{x}_{s}\widetilde{x}_{s}}+\boldsymbol{\Phi}\mathbf{E}_{f_{s}f_{s}}\boldsymbol{\Phi}^{t}\right.\right.\right.\\\
&\left.\left.\left.+\boldsymbol{\Lambda}_{s}\mathbf{E}_{l_{s}l_{s}}\boldsymbol{\Lambda}_{s}^{\top}-2\mathbf{E}_{\widetilde{x}_{s}f_{s}}\boldsymbol{\Phi}^{\top}-2\mathbf{E}_{\widetilde{x}_{s}l_{s}}\boldsymbol{\Lambda}_{s}^{\top}+2\boldsymbol{\Phi}\mathbf{E}_{f_{s}l_{s}}\boldsymbol{\Lambda}_{s}^{\top}\right)\right]\right),\end{split}$
(7)

with
$\mathbf{E}_{\widetilde{x}_{s}\widetilde{x}_{s}},\mathbf{E}_{f_{s}f_{s}},\mathbf{E}_{l_{s}l_{s}},\mathbf{E}_{\widetilde{x}_{s}f_{s}}$
and $\mathbf{E}_{f_{s}l_{s}}$ as in Algorithm 1. We obtain the updates for
each parameter setting the partial derivatives of $l_{c}(\boldsymbol{\theta})$
to zero, and update each of them iteratively while keeping the others fixed at
their most recent update. The updates for the precision
$\boldsymbol{\Psi}_{s}$, common $\boldsymbol{\Phi}$, and study-specific
loadings matrices $\boldsymbol{\Lambda}_{s}$ are obtained analogously to the
ECM updates in the MSFA (De Vito et al., 2019), with
${\widetilde{\mathbf{x}}_{is}}=\mathbf{x}_{is}{\color[rgb]{0,0,0}-\boldsymbol{\beta}\mathbf{b}_{is}}$,
maximizing Equation (7) w.r.t.
$\boldsymbol{\Phi},\boldsymbol{\Lambda}_{s},\boldsymbol{\Psi}_{s}$:

$\begin{split}\widehat{\boldsymbol{\Psi}}_{s}=&\text{diag}\left\\{\mathbf{E}_{\widetilde{x}_{s}\widetilde{x}_{s}}+\widehat{\boldsymbol{\Phi}}\mathbf{E}_{ff}\widehat{\boldsymbol{\Phi}}^{\top}+\widehat{\boldsymbol{\Lambda}}_{s}\mathbf{E}_{l_{s}l_{s}}\widehat{\boldsymbol{\Lambda}}_{s}^{\top}-2\mathbf{E}_{\widetilde{x}_{s}f}\widehat{\boldsymbol{\Phi}}^{\top}-\Biggl{.}2\mathbf{E}_{\widetilde{x}_{s}l_{s}}\widehat{\boldsymbol{\Lambda}}_{s}^{\top}+2\widehat{\boldsymbol{\Phi}}\mathbf{E}_{fl_{s}}\widehat{\boldsymbol{\Lambda}}_{s}^{\top}\right\\}\\\
\mbox{vec}(\widehat{\boldsymbol{\Phi}})=&\sum_{s=1}^{S}\left(\mathbf{E}_{ff}^{\top}\otimes
n_{s}\widehat{\boldsymbol{\Psi}}^{-1}_{s}\right)\mbox{vec}\left(n_{s}\widehat{\boldsymbol{\Psi}}^{-1}_{s}\mathbf{E}_{\widetilde{x}_{s}f}-n_{s}\widehat{\boldsymbol{\Psi}}^{-1}_{s}\widehat{\boldsymbol{\Lambda}}_{s}\mathbf{E}_{fl_{s}}^{\top}\right)\\\
\widehat{\boldsymbol{\Lambda}}_{s}=&\Biggl{(}\mathbf{E}_{x_{s}f}-\widehat{\boldsymbol{\Phi}}\mathbf{E}_{fl_{s}}\Biggr{)}\Biggl{(}\mathbf{E}_{l_{s}l_{s}}\Biggr{)}^{-1},\end{split}$
(8)

where $\otimes$ is the Kronecker product, vec is the vec operator, and the
linear equation is solved with the Lyapunov Equation. More details about the
optimizations can be found in Supplementary A.

To estimate the coefficients $\boldsymbol{\beta}$, we first expand Equation
(3):

$\begin{split}l_{c}(\boldsymbol{\theta})\propto&\sum_{s=1}^{S}\Biggl{\\{}-\frac{n_{s}}{2}tr\left[\boldsymbol{\Psi}_{s}^{-1}\frac{\sum_{i}^{n_{s}}\left(-2\mathbf{x}_{is}\mathbf{b}^{\top}_{is}\boldsymbol{\beta}^{\top}+\boldsymbol{\beta}\mathbf{b}_{is}\mathbf{b}_{is}^{\top}\boldsymbol{\beta}^{\top}+2\boldsymbol{\beta}\mathbf{b}_{is}\mathbf{E}_{f_{i}}^{\top}\boldsymbol{\Phi}^{\top}+2\boldsymbol{\beta}\mathbf{b}_{is}\mathbf{E}_{l_{is}}^{\top}\boldsymbol{\Lambda}_{s}^{\top}\right)}{n_{s}}\right]\Biggr{\\}}.\\\
\end{split}$ (9)

We then take the derivative and the expected values of Equation (9)

$\begin{split}\frac{\partial
l_{c}(\boldsymbol{\theta})}{\partial\boldsymbol{\beta}}=&\sum_{s=1}^{S}\left(-\frac{n_{s}}{2}\boldsymbol{\Psi}_{s}^{-1}\left[\frac{\sum_{i}^{n_{s}}2\left(\boldsymbol{\beta}\mathbf{b}_{is}\mathbf{b}_{is}^{\top}-(\mathbf{x}_{is}-\boldsymbol{\Phi}\mathbf{E}_{f_{i}}-\boldsymbol{\Lambda}_{s}\mathbf{E}_{l_{is}}\right)\mathbf{b}_{is}^{\top}}{n_{s}}\right]\right)=0,\end{split}$
(10)

and solve Equation (10) to obtain:

$\begin{split}\widehat{\boldsymbol{\beta}}=&\left[\sum_{s=1}^{S}\sum_{i}^{n_{s}}(\mathbf{x}_{is}-\widehat{\boldsymbol{\Phi}}\mathbf{E}_{f_{i}}-\widehat{\boldsymbol{\Lambda}}_{s}\mathbf{E}_{l_{is}})\mathbf{b}_{is}^{\top}\right]\left[\sum_{s=1}^{S}\sum_{i}^{n_{s}}(\mathbf{b}_{is}\mathbf{b}_{is}^{\top})\right]^{-1}.\end{split}$
(11)

The stopping criterion used for our proposed ECM algorithm is based on Aitken
acceleration (McLachlan, 2008). Algorithm 1 is stopped when a tolerance
$\epsilon^{*}=\mid l^{(t+1)}-l^{(t)}\mid$, with
$l^{(t+1)}=l_{c}(\boldsymbol{\theta}^{(t)})+\frac{l_{c}(\boldsymbol{\theta}^{(t+1)})-l_{c}(\boldsymbol{\theta}^{(t)})}{1-c^{(t)}}$,
and
$c^{(t)}=\frac{l_{c}(\boldsymbol{\theta}^{(t+1)})-l_{c}(\boldsymbol{\theta}^{(t)})}{l_{c}(\boldsymbol{\theta}^{(t)})-l_{c}(\boldsymbol{\theta}^{(t-1)})}$
in the complete log-likelihood change is reached, or a maximum number of
iterations $T$ have been completed. By default we set $\epsilon^{*}=10^{-7}$
and $T=50,000$.

The selection of the dimension of the common and study-specific factors $q$
and $q_{s}$ is made by choosing the model that leads to the best Akaike
information criterion (AIC) (Akaike, 1974) or the Bayesian information
criterion (BIC) (Schwarz, 1978). More details about the latent cardinality
selection will be discussed in Sections §4 and §5.

### Initialization of the ECM algorithm

Parameter initialization is crucial to obtain better local modes and reduce
computational time, as the ECM algorithm can be sensitive to parameter
initialization. Here, we propose a simple two-step least-squares
initialisation by stacking all the $S$ study-data in a single data set
$\mathbf{X}^{\top}\in\mathbb{R}^{n\times p},n=n_{1}+\dots+n_{S}$:

  * Step 1

initialize
$\boldsymbol{\beta}^{(0)}=\left[\mathbf{B}\mathbf{B}^{\top}\right]^{-1}\mathbf{B}\mathbf{X}^{\top}$
with $\mathbf{X}^{\top}\in\mathbb{R}^{n\times
p},\mathbf{B}^{\top}\in\mathbb{R}^{n\times p_{b}}$

  * Step 2

let $\widetilde{\mathbf{X}}=\mathbf{X}-\boldsymbol{\beta}^{(0)}\mathbf{B}$.
Initialize the extra parameters as the Multi-study factor analysis:

    1. 1.

perform a Principal Components Analysis (PCA) on $\widetilde{\mathbf{X}}$ and
initialized $\boldsymbol{\Phi}^{(0)}$ with the first $q$ principal components,
and

    2. 2.

perform factor analysis for each $\widetilde{\mathbf{X}}_{s},s=1,\dots,S$, and
initialize $\boldsymbol{\Lambda}^{(0)}$ and $\boldsymbol{\Psi}^{(0)}$ with the
obtained loading matrices and unharnesses, respectively.

### Identifiability Considerations

The MSFR model, like other standard factor analysis models, is non-
identifiable up orthogonal transformations of the loadings (common
$\boldsymbol{\Phi}$ and study specific $\boldsymbol{\Lambda}_{s}$) and
loadings sign switch. In addition to those sources of non-identification, the
MSFR model has an indeterminacy of the covariance $\boldsymbol{\Sigma}_{s}$
decomposition.

The rotation identifiability arises since any orthogonal transformation of the
factors and the loading matrices, common and study-specific, lead to the same
distribution of $\mathbf{x}_{is}$. Thus, the marginal distribution of
$\mathbf{x}_{is}$ can equivalently be redefined as a multivariate normal with
covariance matrix
$\boldsymbol{\Sigma}_{s}=\widetilde{\boldsymbol{\Phi}}\widetilde{\boldsymbol{\Phi}}^{\top}+\widetilde{\boldsymbol{\Lambda}}_{s}\widetilde{\boldsymbol{\Lambda}}_{s}^{\top}+\boldsymbol{\Psi}_{s}$,
where $\widetilde{\boldsymbol{\Phi}}=\boldsymbol{\Phi}\mathbf{R}$ and
$\widetilde{\boldsymbol{\Lambda}}_{s}=\boldsymbol{\Lambda}_{s}\mathbf{R}_{s}$,
with $\mathbf{R}\in\mathbb{R}^{q\times q}$ and
$\mathbf{R}_{s}\in\mathbb{R}^{q_{s}\times q_{s}}$ orthogonal matrices.
Furthermore, Model (2) can be rewritten as
$\mathbf{x}_{is}={\boldsymbol{\beta}}\mathbf{b}_{is}+\widetilde{\boldsymbol{\Phi}}\widetilde{\mathbf{f}}_{is}+\widetilde{\boldsymbol{\Lambda}}_{s}\widetilde{\mathbf{l}}_{is}+\mathbf{e}_{is}$,
with $\widetilde{\mathbf{f}}_{is}=\mathbf{R}^{\top}{\mathbf{f}}_{is}$ and
$\widetilde{\mathbf{l}}_{is}=\mathbf{R}_{s}^{\top}\mathbf{l}_{is}$. In the
factor analysis literature, several strategies have been developed to tackle
this issue, such as performing a varimax rotation (Kaiser, 1958), or
restricting the loading matrix to be lower triangular (Lopes and West, 2004;
Carvalho et al., 2008) or generalized lower triangular (Frühwirth-Schnatter et
al., 2023). Here, we extend the traditional varimax rotation to both the
$\boldsymbol{\Phi}$ and $\boldsymbol{\Lambda}_{s}$. Such alternative was
selected due to its interpretability. MSFR is also not uniquely identified due
to sign switching of the loading matrices $\boldsymbol{\Phi}$ and
$\boldsymbol{\Lambda}_{s}$. We fix this mild issue by constraining the sign of
a subset of loadings.

Under varimax rotation and sign constraints, a full rank loading matrix is
uniquely determined in regular factor analysis, giving the error term.
However, in multi-study latent models, an additional source of non-
identifiability arises, the indeterminacy of the covariance
$\boldsymbol{\Sigma}_{s}$ decomposition. More in detail, the difference
between the covariance matrix and the idiosyncratic precisions can be
redefined as
$\boldsymbol{\Sigma}_{s}-\boldsymbol{\Psi}_{s}={\boldsymbol{\Phi}^{*}}{\boldsymbol{\Phi}}^{*\top}+{\boldsymbol{\Lambda}}^{*}_{s}{\boldsymbol{\Lambda}}_{s}^{*\top}$,
from where column vectors of $\boldsymbol{\Phi}$ and
${\boldsymbol{\Lambda}}_{s}$ can move from one loading matrix to the other. To
solve this issue, we follow the strategy of (De Vito et al., 2019), and
constraint the
$[\boldsymbol{\Phi},\boldsymbol{\Lambda}_{1},\dots,\boldsymbol{\Lambda}_{S}]\in\mathbb{R}^{p\times(q+\sum_{s=1}^{S}q_{s})}$
matrix, containing both the common and study specific loadings, to have a full
column rank $q+\sum_{s=1}^{S}q_{s}\leq p$. Thus, the span of the column
vectors of $\boldsymbol{\Sigma}_{s}-\boldsymbol{\Psi}_{s}$ is equal to the sum
of the span of the column vectors of $\boldsymbol{\Phi}$ and
${\boldsymbol{\Lambda}}_{s}$, leading to uniquely determined loading matrices.

##  4 Simulation Study

We first use simulation studies to assess the ability of our MSFR model to
detect the true signals for both the common and the study-specific and
outperform the standard MSFA in the presence of observed covariates.
Specifically, we compare the MSFR with MSFA to illustrate how much inference
can suffer when not accounting for covariate effects. To evaluate if the
covariate effect can be removed from the data and then adopt the MSFA to the
residuals, we fit a linear regression (LR) model, followed by MSFA estimation
on the residuals, here called (MSFA&LR). Moreover, we compare the MSFR with
the vanilla Factor Regression (FR) (Avalos-Pacheco et al., 2022; Carvalho et
al., 2008) to assess the advantages of estimating the study-specific factors.
The R-code for our methods is available at https://github.com/rdevito/MSFR.

We consider various scenarios with different sample sizes $n_{s}$, number of
studies $S$, number of observed variables $p$, structure of
$\boldsymbol{\beta}$, and structure of $\boldsymbol{\Sigma}$, as shown in
Table 2. Scenario 1 shows the advantages of our method in situations with
small covariate effects. Scenarios 2 and 3 closely mimic the case study data,
the HCHS/SOL dataset, with its same number of studies and patient
characteristics varying only in the structure of $\boldsymbol{\beta}$,
changing the dimension of $p_{b}$. For all the scenarios, we simulated 100
collections of datasets from
$N(\boldsymbol{\beta}\mathbf{b}_{is},\boldsymbol{\Sigma}_{s})$, with
$\boldsymbol{\Sigma}_{s}=(\boldsymbol{\Phi}\boldsymbol{\Phi}^{\top}+\boldsymbol{\Lambda}_{s}\boldsymbol{\Lambda}_{s}^{\top}+\boldsymbol{\Psi}_{s})$.
The common $\boldsymbol{\Phi}$ and study specific $\boldsymbol{\Lambda}_{s}$
loading matrices are assumed to be sparse, with one-third of non-zero entries,
randomly generated from the uniform densities
$\boldsymbol{\Phi}\sim\text{sign}\times U(0.6,1)$, with $\text{sign}\sim
U\\{-1,1\\}$ and $\boldsymbol{\Lambda}_{s}\sim U(-1,1)$. Finally, the
idiosyncratic precision matrices, $\boldsymbol{\Psi}_{s}$, were assumed to be
diagonal matrices with the entries generated uniformly between (0,1).

Table 2 shows the mean results across 100 simulations of (i) the selected
number of common and study-specific factors, $\widehat{q}$ and
$\widehat{q}_{s}$ respectively, and (ii) the RV coefficient (Robert and
Escoufier, 1976) between the ground truth and the reconstructed value for the
coefficients $\widehat{\boldsymbol{\beta}}$, the common
$\widehat{\boldsymbol{\Phi}}$ and the study-specific loadings
$\widehat{\boldsymbol{\lambda}}_{s}$, and the total covariance
$\widehat{\boldsymbol{\Sigma}}_{s}$ for all the studies. The RV coefficient
$\in[0,1]$ is a multivariate generalization of the squared Pearson correlation
coefficient. The higher the value of the RV coefficient, the better the
reconstruction. For each simulation, we considered a grid of possible values
for $\widehat{q}$ and $\widehat{q}_{s}$, and we fit the models for all the
grid values and performed model selection. Table 2 displays the results for
the selected models following the Akaike information criterion (AIC) (Akaike,
1974) and the Bayesian information criterion (BIC) (Schwarz, 1978) approaches
for the three different scenarios.

The results showed the advantages of using our proposed method. MSFR presented
the highest RV for the common and study-specific loading matrices as well as
for the whole covariances. We can also see that MSFR estimated $\widehat{q}$
and $\widehat{q}_{s}$ accurately in all the scenarios, independently of the
model selection criteria (BIC and AIC), showing a most robust behaviour than
MSFA&LR and MSFA. As expected, FR required a higher number of factors to take
into account the impact of the study-specific factors. Likewise, MSFA required
more factors since it does not consider the covariate effects. MSFA&LR
outperformed all the methods in the reconstruction of $\boldsymbol{\beta}$,
but showed the worst performance in the reconstruction of $\widehat{q}$,
$\widehat{q}_{s}$, $\boldsymbol{\Lambda}_{s}$ and $\boldsymbol{\Phi}$ in all
the scenarios, highlighting the limitations of using two-step procedures
relative to a joint estimation of the multiple-factor model and covariate
effects. We then assess the reconstructions of $\boldsymbol{\beta}$,
$\boldsymbol{\Sigma}_{\boldsymbol{\Phi}}$, $\boldsymbol{\lambda}_{s}$ and
$\boldsymbol{\Sigma}_{s}$ obtained with the AIC criteria for all the methods
and scenarios.

Table 2: Average RV over 100 simulations

|  | $\widehat{q}$ | $\widehat{q}_{s}$ | $\widehat{\boldsymbol{\beta}}$ | $\widehat{\boldsymbol{\Phi}}$ | $\widehat{\boldsymbol{\lambda}}_{1}$ | $\widehat{\boldsymbol{\lambda}}_{2}$ | $\widehat{\boldsymbol{\lambda}}_{3}$ | $\widehat{\boldsymbol{\lambda}}_{4}$ | $\widehat{\boldsymbol{\lambda}}_{5}$ | $\widehat{\boldsymbol{\lambda}}_{6}$ | $\widehat{\boldsymbol{\Sigma}}_{1}$ | $\widehat{\boldsymbol{\Sigma}}_{2}$ | $\widehat{\boldsymbol{\Sigma}}_{3}$ | $\widehat{\boldsymbol{\Sigma}}_{4}$ | $\widehat{\boldsymbol{\Sigma}}_{5}$ | $\widehat{\boldsymbol{\Sigma}}_{6}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Scenario 1: $q=3,q_{s}=1,S=2,p_{b}=2,p=20,n_{s}=\\{500,500,500\\}$
BIC | MSFR | 3.00 | 1.00 | 0.978 | 0.996 | 0.959 | 0.931 |  |  |  |  | 0.974 | 0.987 |  |  |  |
MSFA&LR | 3.00 | 1.00 | 0.993 | 0.973 | 0.922 | 0.875 |  |  |  |  | 0.923 | 0.942 |  |  |  |
MSFA | 4.95 | 1.04 | NA | 0.915 | 0.900 | 0.844 |  |  |  |  | 0.889 | 0.904 |  |  |  |
FR | 3.57 | NA | 0.991 | 0.956 | NA | NA |  |  |  |  | 0.912 | 0.888 |  |  |  |
AIC | MSFR | 3.14 | 1.19 | 0.978 | 0.994 | 0.936 | 0.904 |  |  |  |  | 0.972 | 0.986 |  |  |  |
MSFA&LR | 1.51 | 3.20 | 0.993 | 0.626 | 0.365 | 0.178 |  |  |  |  | 0.585 | 0.608 |  |  |  |
MSFA | 2.80 | 3.94 | NA | 0.759 | 0.399 | 0.188 |  |  |  |  | 0.719 | 0.736 |  |  |  |
FR | 3.84 | NA | 0.991 | 0.948 | NA | NA |  |  |  |  | 0.913 | 0.885 |  |  |  |
Scenario 2:
$q=4,q_{s}=1,S=6,p_{b}=7,p=42,n_{s}=\\{1257,1444,2126,4940,2314,897\\}$
BIC | MSFR | 4.00 | 1.00 | 0.952 | 1.000 | 0.983 | 0.982 | 0.990 | 0.995 | 0.991 | 0.981 | 0.999 | 0.999 | 0.999 | 0.999 | 0.999 | 0.999
MSFA&LR | 1.14 | 4.13 | 0.999 | 0.602 | 0.017 | 0.036 | 0.033 | 0.002 | 0.006 | 0.011 | 0.841 | 0.848 | 0.858 | 0.866 | 0.832 | 0.847
MSFA | 8.22 | 2.51 | NA | 0.953 | 0.927 | 0.908 | 0.939 | 0.982 | 0.973 | 0.927 | 0.965 | 0.965 | 0.964 | 0.966 | 0.966 | 0.965
FR | 5.00 | NA | 0.998 | 0.971 | NA | NA | NA | NA | NA | NA | 0.860 | 0.874 | 0.868 | 0.872 | 0.876 | 0.892
AIC | MSFR | 4.12 | 1.00 | 0.952 | 1.000 | 0.983 | 0.982 | 0.990 | 0.995 | 0.991 | 0.981 | 0.999 | 0.999 | 0.999 | 0.999 | 0.999 | 0.999
MSFA&LR | 3.00 | 5.00 | 0.999 | 0.815 | 0.206 | 0.235 | 0.222 | 0.207 | 0.208 | 0.218 | 0.931 | 0.932 | 0.937 | 0.938 | 0.930 | 0.932
MSFA | 8.39 | 5.00 | NA | 0.953 | 0.924 | 0.907 | 0.938 | 0.981 | 0.971 | 0.929 | 0.966 | 0.966 | 0.965 | 0.967 | 0.966 | 0.965
FR | 5.00 | NA | 0.998 | 0.971 | NA | NA | NA | NA | NA | NA | 0.860 | 0.874 | 0.868 | 0.872 | 0.876 | 0.892
Scenario 3:
$q=4,q_{s}=1,S=6,p_{b}=9p=42,n_{s}=\\{1257,1444,2126,4940,2314,897\\}$
BIC | MSFR | 3.99 | 1.00 | 0.940 | 0.999 | 0.970 | 0.975 | 0.982 | 0.986 | 0.981 | 0.971 | 0.999 | 0.999 | 0.999 | 0.999 | 0.999 | 0.999
MSFA&LR | 1.36 | 4.17 | 0.998 | 0.617 | 0.033 | 0.074 | 0.050 | 0.022 | 0.026 | 0.030 | 0.852 | 0.856 | 0.864 | 0.869 | 0.847 | 0.851
MSFA | 8.06 | 2.28 | NA | 0.923 | 0.925 | 0.903 | 0.936 | 0.982 | 0.976 | 0.920 | 0.929 | 0.930 | 0.929 | 0.929 | 0.931 | 0.928
FR | 5.00 | NA | 0.998 | 0.971 | NA | NA | NA | NA | NA | NA | 0.859 | 0.874 | 0.867 | 0.871 | 0.876 | 0.892
AIC | MSFR | 4.16 | 1.00 | 0.940 | 0.999 | 0.970 | 0.975 | 0.982 | 0.986 | 0.981 | 0.971 | 0.999 | 0.999 | 0.999 | 0.999 | 0.999 | 0.999
MSFA&LR | 2.60 | 5.00 | 0.998 | 0.787 | 0.126 | 0.181 | 0.136 | 0.120 | 0.120 | 0.135 | 0.923 | 0.923 | 0.929 | 0.928 | 0.923 | 0.922
MSFA | 9.06 | 5.00 | NA | 0.923 | 0.923 | 0.898 | 0.936 | 0.982 | 0.975 | 0.922 | 0.928 | 0.929 | 0.928 | 0.928 | 0.930 | 0.927
FR | 5.00 | NA | 0.998 | 0.971 | NA | NA | NA | NA | NA | NA | 0.859 | 0.874 | 0.867 | 0.871 | 0.876 | 0.892

Supplementary B Figure S1 displays the heatmaps of the true
$\boldsymbol{\beta}$ and $\boldsymbol{\Sigma}_{\boldsymbol{\Phi}}$, and their
average estimations (left panel), as well as the RV’s box-plots for
$\boldsymbol{\beta}$ and ${\boldsymbol{\Phi}}$ (right panel). Figure S2 shows
the box-plots of the RV coefficients of $\boldsymbol{\lambda}_{s}$ and
$\boldsymbol{\Sigma}_{s}$. Overall, MSFR performs better than the other
methods, with the smallest variance in all cases, except for
$\boldsymbol{\beta}$ (Figure S1). It is important to highlight that even
though MSFA achieved a competitive reconstruction of $\boldsymbol{\Phi}$, the
average $\boldsymbol{\Sigma}_{\boldsymbol{\Phi}}$ estimates over the
simulations are less precise, as displayed in the heatmaps in Figure S2.
Finally, Supplementary C, Figures S3 and S4 present the results of the
reconstructions of $\boldsymbol{\beta}$,
$\boldsymbol{\Sigma}_{\boldsymbol{\Phi}}$, $\boldsymbol{\lambda}_{s}$ and
$\boldsymbol{\Sigma}_{s}$ obtained via the BIC criterion. The results showed
that the performance of our competitors varied significantly, further
corroborating that those methods may not be robust to the model selection
criteria.

##  5 The HCHS/SOL application

We demonstrate the usefulness of our method by performing four different
analyses on the HCHS/SOL. First, we compare the prediction accuracy between
the MSFR and the MSFA to show the power of our method in predictive ability.
Second, we explore the common and study-specific factors loadings among the
five different background groups of the HCHS/SOL. Third, we associate the
common and study-specific loadings with food groups. Fourth, we validate our
obtained dietary patterns, both common and study-specific, associating each
factor with the alternative eating index(AHEI-2010) - a measure of overall
diet quality previously related to cardiometabolic disease risk (Liese et al.,
2015).

Predictive analyses. We compare the k-fold cross-validation (setting k=5)
prediction errors computed by our MSFR and the MSFA. Thus, we fit both the
MSFR and the MSFA model on a random 80% of the data, followed by a different
80% for different folds. We then evaluate the prediction error on the
remaining 20% in each fold. Predictions are performed as

$\mbox{MSFR:
}\widehat{\mathbf{x}}_{is}=\widehat{\boldsymbol{\beta}}\widehat{\mathbf{b}}_{is}+\widehat{\boldsymbol{\Phi}}^{MSFR}\widehat{\mathbf{f}}_{is}+\widehat{\boldsymbol{\Lambda}}^{MSFR}_{s}\widehat{\mathbf{l}}_{is}\;\;\;\;\;\;\;\;\;\;\;\;\mbox{
MSFA:
}\widehat{\mathbf{x}}_{is}=\widehat{\boldsymbol{\Phi}}^{MSFA}\widehat{\mathbf{f}}_{is}+\widehat{\boldsymbol{\Lambda}}^{MSFA}_{s}\widehat{\mathbf{l}}_{is},$

where $\widehat{\boldsymbol{\Phi}}^{MSFR}$ is the common factor loadings
obtained with the MSFR method, while $\widehat{\boldsymbol{\Phi}}^{MSFA}$ is
the common factor loadings obtained with the MSFA method;
$\widehat{\boldsymbol{\Lambda}}^{MSFR}_{s}$ is specific factor loadings
estimated with MSFR, while $\widehat{\boldsymbol{\Lambda}}^{MSFA}_{s}$ is
specific factor loadings. The factor scores are then computed using the
Bartlett method since scores are uncorrelated between each other (DiStefano et
al., 2009), see Supplement D for further details, by using different
parameters obtained in the MSFR and MSFA. Thurstone method for estimating
factor score lead to similar results. We then assessed the mean squared error
(MSE) of prediction for each fold as
$MSE=\frac{1}{n}\sum_{s=1}^{S}\sum_{i=1}^{n_{s}}(\mathbf{x}_{is}-\widehat{\mathbf{x}}_{is})^{2},$
using the 20% of the samples in each study set aside for each cross-validation
iteration. We averaged all the MSE obtained in each fold. Table 3 shows the
obtained MSE for the MSFR and the MSFA, obtained with the factor scores
computed via Bartlett and Thurstone. The MSE is 4% smaller for MSFR than for
MSFA. This analysis reveals how MSFR borrows strength across studies in
estimating the factor loadings and clearly detecting factor loading signals by
removing all the covariates noise that can impact them. Thus, our predictive
power in independent observations is not only preserved but also improved.

Method | MSE
---|---
Bartlett | Thurstone
MSFR | 1.70 | 1.62
MSFA | 1.76 | 1.68

Table 3: Mean square error (MSE) obtained via k-fold cross validation for both
the multi-study factor regression model (MSFR) and the multi-study factor
analysis model (MSFA). Factor scores are obtained with the Bartlett and
Thurstone method.

Factor loadings analyses.

We considered two main tasks: to analyze the estimated common and study-
specific factors loadings, $\boldsymbol{\Phi}$ and $\boldsymbol{\Lambda}_{s}$,
respectively, and to study the loading association with different food groups.

The MSFR model selected $\widehat{q}=4$ common factors and $\widehat{q}_{s}=1$
study-specific factor for the six different backgrounds (obtained with the AIC
procedure). The common factors explained 62.5% of the total variation, each of
which explained 23.8%, 15.5%, 9.8% and 6.9% of the total variance. The total
variance explained by each study-specific factor range between 6.9%
(Dominican, Mexican, and Puerto Rican background) to 6.5% (Cuban background).
Thus, no differences are detected among the explained variance for each study-
specific loading. We remark that all the factors, common and study-specific,
explain more than 5% of the total variance, a trade-off generally considered
in factor analysis, confirming our factor dimension estimation.

We proceed by inspecting the factor structures. Figure 1 provides a visual
representation of the common and study-specific loadings, on the left and
right panels, respectively. In order to facilitate interpretation and thus
name each factor, we performed a varimax rotation (Kaiser, 1958) for the
common loading matrix. This rotation was unnecessary for
$\boldsymbol{\Lambda}_{s}$ as our MSFR estimate only one extra study-specific
loading.

The varimax rotation in the common loadings maximizes the variance in the
squared rotated loadings, making the loadings extreme and positive for a
better interpretation. Factor 1, namely Plant-based foods, showed high-valued
loadings on vegetable protein, phosphorus, magnesium, iron, zinc, copper,
potassium, manganese, thiamin, niacin, pantothenic acid, vitamin B6, natural
folate, soluble and insoluble dietary fiber. Factor 2, namely Animal products,
displayed high-valued loadings in animal protein, cholesterol, long-chain
saturated fatty acids (LCSFA), long-chain monounsaturated fatty acids (LCMFA),
linolenic acid, arachidonic acid, total trans fatty acid, selenium, sodium,
and natural alpha-tocopherol. Factor 3, namely Dairy products, exhibited high-
valued loadings on short- and medium-chain saturated fatty acids, calcium,
riboflavin, retinol, and vitamin D. Factor 4, namely Seafood, presented high-
valued loadings on eicosapentaenoic acid (EPA), docosapentaenoic acid (DPA),
and docosahexaenoic acid (DHA).

In the absence of any rotation in the specific factor loadings, we based our
interpretation on the opposite (positive and negative) signs or colors (blue
and dark red, respectively) illustrated in the study-specific factor loading
heatmap (Figure 1 right panel). They displayed some similarities, showing a
general behavior with positive loadings for the animal food source (i.e.,
animal proteins, cholesterol, arachidonic acid, vitamin B12, vitamin D) and
negative loadings for the vegetable source (i.e., vegetable protein, starch,
and natural folate). However, the study-specific factor loadings presented
some differences across each other (Figure 1 right panel, in black box),
providing a deeper understanding of the background-specific dietary patterns.
In particular, the short-, medium-, long-chain saturated fatty acids and the
long-chain monounsaturated fatty acids are more expressed for the Central
America pattern than the other background patterns, revealing that the Central
America background is more associated with the good oily components. In fact,
these nutrients can be found in nuts, seeds, and coconut oil. The Cuban
background pattern expressed higher loadings in calcium, phosphorus, and
thiamin than the other background groups but not in potassium. This suggests
that the Cuban background pattern is enriched with more vegetable source foods
than the other background groups. Finally, the Dominican and Mexican
background patterns showed higher-valued loadings for the fiber (soluble and
insoluble) and not in the pantothenic acid and vitamin B6. This suggests that
the Dominican and Mexican background patterns can be associated with grain
consumption. The other two background patterns, i.e., Puerto-Rican and South
American, are more similar, representing some variants of an animal profile.

Figure 1: Heatmap of the common (left panel) and study-specific (right panel)
factor loadings in the HCHS/SOL.

Associations with food groups.

To further corroborate and validate our results, we associate each common and
study-specific factor with the food groups from the 24-hour recalls. Thus, we
estimate the factor scores, i.e., the degree to which each subject’s diet
adheres to one of the identified factors. In order to do that, we adopted the
Bartlett and Thurstone approaches (Johnson, 1988; DiStefano et al., 2009) in
the MSFR model. The correlations between Bartlett and Thurstone scores
exceeded 0.95 for the common and study-specific factors, so we adopted the
Bartlett method since scores are uncorrelated between each other (DiStefano et
al., 2009), see Supplement D for further details. We reported the association
between the factors, common and study-specific, and the food group in Figure
2, left and right panels, respectively, where in grey is reported the high
consumption of the corresponding food in that factor.

Figure 2: Association of each factor with each food group for each common
factor (first box) and each study-specific factor (second box).

The names attributable to each common factor, i.e., Plant-based food, Animal
products, Dairy products, and Seafood, were further confirmed by this
analysis. In fact, the Plant-based foods pattern is highly enriched in all
types of fruits (Fruit-Citrus, Fruit-Other, Fruit-All) and vegetarian (Dark
Green and Veg-Tomatoes). This pattern is also enriched in grain, especially
refined and soy, potentially related to the Hispanic diet, with corn, soybean,
and tortilla. The Animal products pattern is enriched in Deli-Meat, Organ
Meats, Meat-All, Cheese, Dairy Dessert and Milk. The pattern of Dairy products
is enriched in all the dairy foods, i.e., Milk, Cheese, Yogurt, Dairy Dessert,
and Milk. These last two dietary patterns express a crucial distinction from
one another, one more associated with animal fat and animal consumption, the
other more oriented on dairy consumption. The Seafood dietary pattern is
enriched in Fish but also in Cheese, Dairy dessert. It is the only pattern to
be enriched in alcohol consumption.

The study-specific analyses confirm our first visual inspection. All the
background patterns reported high meat, red, poultry, and organ meats
consumption. Two background patterns, i.e., Puerto-Rican and South American,
are similar, with high consumption in all the meat types (including poultry),
with some nuances. For example, the Puerto-Rican pattern presents high
consumption of alcohol compared to the South Americans. We called Puerto-Rican
and South American patterns the “animal source-only” patterns. High
consumption of alcohol is also presented in two other background dietary
patterns, Central American and Cuban. The Central American pattern is very
similar to the “animal source-only” patterns, with some slight differences,
i.e., it does not show the deli-meat consumption and is enriched in milk. The
Dominican and Mexican patterns present high consumption of grain, especially
refined grain and grain-all, revealing these two patterns to be not only
animal sources. The Cuban pattern is the only one to show high consumption of
vegetables (Dark green) and potatoes together with whole grains. Thus, this
pattern is based on vegetable sources.

Association with the Alternative Healthy Eating Index(AHEI-2010).

We further performed an analysis to validate our dietary patterns, associating
both the common and the study-specific patterns with the Alternative Healthy
Eating Index (AHEI-2010). The AHEI-2010 is a metric to determine if a group of
foods, and so nutrients, adheres to the Dietary Guidelines for Americans
(DGA)(Chiuve et al., 2012). The alternative eating index was adjusted based on
foods and nutrients preventing chronic disease risk. The AHEI-2010 is an
11-component diet quality score ranging from 0 to 110. The higher the
AHEI-2010 index, the lower the chronic disease risk. Thus, it is a crucial
measure to detect if patterns can be related to any possible risk, especially
cardiometabolic disease risk (Liese et al., 2022).

We assessed the association between our dietary patterns and the AHEI-2010 by
estimating the mean score within quantile-based categories of factor scores
using a weighted generalized linear regression model adjusted for total energy
intake only. We computed the mean and the standard error for each estimate
(Figure 3). We remark that we have not included any covariates as confounders
since we have already accounted for those in our multi-study factor regression
model.

For the common dietary patterns, we observed that both the Plant-based foods
and the Seafood patterns increased the AHEI-2010, while both the Animal
Products and Dairy products patterns decreased the AHEI-2010. The Animal
Products pattern is strongly decreasing the AHEI-2010, scoring from 49.7 (
first quantile) to 45.2 ( fifth quantile), while the Dairy products pattern is
slightly decreasing, scoring from 48.1 ( first quantile) to 46.7 (fifth
quantile).

The study-specific patterns have similar behavior, all decreasing the AHEI,
with some differences. For instance, the Dominican background dietary pattern
increases strongly the AHEI, ranging from 51.2 (first quantile) to 46.4 (third
quantile). The Cuban background dietary pattern associated with more
vegetables slightly decreases the AHEI, ranging from 44.7 (first quantile) to
43.5 (third quantile). The standard errors for the study-specific are all very
low, close to zero.

Figure 3: Mean (top panels) association of each factor with the AHEI-2010 and
its standard deviation (bottom panels), for both the common (left side of the
dashed line) and the study-specific (right side of the dashed line)

##  6 Discussion

In this article, we aim to identify the shared and ethnic background site-
specific dietary patterns in the HCHS/SOL study, while adjusting for socio-
demographic characteristics or confounders that do not affect the dietary
patterns like: gender, education level and alcohol and cigarette intake and
BMI. To achieve this goal, we have presented a novel class of factor models:
the Multy-Study Factor Regression (MSFR).

MSFR integrates data from multiple sources using joint dimension reduction,
via Multi-study Factor Analyses (De Vito et al., 2019, 2021), and covariate
effect adjustment, via latent factor regression (Avalos-Pacheco et al., 2022;
Carvalho et al., 2008). Our MSFR is fitted in a deterministic fashion via an
efficient ECM algorithm with closed-form updates and publicly available at
https://github.com/rdevito/MSFR.

Our proposed MSFR has the potential to be adopted in several different data
settings where the goal is to disentangle the similarities and dissimilarities
among heterogeneous populations, studies or groups, with the added benefit of
adjusting for the effects of covariates. The proposed covariate adjustment
essentially makes use of information about relationships between socio-
demographic characteristics and the outcome so that our MSFR can better
identify the common and study-specific loading matrices and latent factors.
Furthermore, MSFR had the most robust estimation of factor cardinality (common
and study-specific), while our competitors MSFA and FR required a higher
number of factors to take into account the covariate effects and the lack of
study specific factors, respectively. The unadjusted MSFA provides a
competitive estimation of the common and study-specific loadings, in cases
when the covariates are weakly prognostic of the outcomes (Scenario 1 of our
simulations). However, their inference may suffer in cases with strongly
prognostic outcomes (Scenarios 2 and 3 of our simulations). Thus, we recommend
the use of our MSFR model generally, regardless of whether our data have a
high or low covariate effect.

In our nutritional epidemiology application, we show that the covariate
adjustment of our MSFR tends to improve the quality of our estimations and
predictions. Further, our model was able to provide interpretable dietary
patterns for the whole Hispanic Community, as well as for each Latino ethnic
background. Our obtained factors, common and study-specific, were explained
and named according to their nutrient values and their correlation to
different food groups. Finally, we illustrated the utility of our method in a
supervised framework by studying the factor association with the alternative
healthy eating index. We emphasize that our proposed method eliminates the
need for an additional covariance adjustment to analyze the AHEI-2010 scores,
which is required in the case of MSFA.

We remark that our novel MSFR can be extended to frameworks of interest beyond
nutritional epidemiology, such as cancer genomics to identify patterns of
common and study-specific transcriptional variation while correcting for pre-
treatment patient characteristics, or finance to study common market trends
and study-specific recession/expansion periods, while adding other economic
variables.

In summary, our model is essential in the case application considered here. We
hope it will help the joint analyses of multiple studies in many different
settings, abating the current challenges of reproducibility of unsupervised
analyses in nutritional epidemiology and data science.

## References

  * Akaike (1974) Akaike, H. (1974).  A new look at the statistical model identification.  IEEE transactions on automatic control 19(6), 716–723.
  * Aune et al. (2016) Aune, D. et al. (2016).  Nut consumption and risk of cardiovascular disease, total cancer, all-cause and cause-specific mortality: a systematic review and dose-response meta-analysis of prospective studies.  BMC medicine 14(1), 1–14.
  * Avalos-Pacheco (2018) Avalos-Pacheco, A. (2018).  Factor regression for dimensionality reduction and data integration techniques with applications to cancer data.  University of Warwick, PhD thesis.
  * Avalos-Pacheco et al. (2022) Avalos-Pacheco, A., D. Rossell, and R. S. Savage (2022).  Heterogeneous large datasets integration using Bayesian factor regression.  Bayesian analysis 17(1), 33–66.
  * Bennett et al. (2022) Bennett, G., L. A. Bardon, and E. R. Gibney (2022).  A comparison of dietary patterns and factors influencing food choice among ethnic groups living in one locality: a systematic review.  Nutrients 14(5), 941.
  * Carvalho et al. (2008) Carvalho, C. M. et al. (2008).  High-dimensional sparse factor modeling: Applications in gene expression genomics.  Journal of the American Statistical Association 103(484), 1438–1456.
  * Castello et al. (2014) Castello, A. et al. (2014).  Spanish mediterranean diet and other dietary patterns and breast cancer risk: case–control epigeicam study.  British journal of cancer 111(7), 1454–1462.
  * Castelló et al. (2022) Castelló, A. et al. (2022).  Adherence to the western, prudent and mediterranean dietary patterns and colorectal cancer risk: Findings from the spanish cohort of the european prospective investigation into cancer and nutrition (epic-spain).  Nutrients 14(15), 3085.
  * Chiuve et al. (2012) Chiuve, S. E. et al. (2012).  Alternative dietary indices both strongly predict risk of chronic disease.  The Journal of nutrition 142(6), 1009–1018.
  * Conway et al. (2009) Conway, D. I. et al. (2009).  Enhancing epidemiologic research on head and neck cancer: Inhance-the international head and neck cancer epidemiology consortium.  Oral oncology 45(9), 743–746.
  * Daviglus et al. (2014) Daviglus, M. L., A. Pirzada, and G. A. Talavera (2014).  Cardiovascular disease risk factors in the hispanic/latino population: lessons from the hispanic community health study/study of latinos (hchs/sol).  Progress in cardiovascular diseases 57(3), 230–236.
  * De Vito et al. (2019) De Vito, R. et al. (2019).  Shared and study-specific dietary patterns and head and neck cancer risk in an international consortium.  Epidemiology (Cambridge, Mass.) 30(1), 93.
  * De Vito et al. (2022) De Vito, R. et al. (2022).  Shared and ethnic background site-specific dietary patterns in the hispanic community health study/study of latinos (hchs/sol).  medRxiv, 2022–06.
  * De Vito et al. (2019) De Vito, R., R. Bellio, L. Trippa, and G. Parmigiani (2019).  Multi-study factor analysis.  Biometrics 75(1), 337–346.
  * De Vito et al. (2021) De Vito, R., R. Bellio, L. Trippa, and G. Parmigiani (2021).  Bayesian multistudy factor analysis for high-throughput biological data.  The annals of applied statistics 15(4), 1723–1741.
  * DiStefano et al. (2009) DiStefano, C., M. Zhu, and D. Mindrila (2009).  Understanding and using factor scores: Considerations for the applied researcher.  Practical Assessment, Research, and Evaluation 14(1), 20\.
  * Edefonti et al. (2012) Edefonti, V. et al. (2012).  Nutrient-based dietary patterns and the risk of head and neck cancer: a pooled analysis in the international head and neck cancer epidemiology consortium.  Annals of Oncology 23(7), 1869–1880.
  * Foster et al. (2003) Foster, G. D. et al. (2003).  A randomized trial of a low-carbohydrate diet for obesity.  New England Journal of Medicine 348(21), 2082–2090.
  * Frühwirth-Schnatter et al. (2023) Frühwirth-Schnatter, S., D. Hosszejni, and H. F. Lopes (2023).  Sparse bayesian factor analysis when the number of factors is unknown.
  * GBD 2013 Risk Factors Collaborators and others (2015) GBD 2013 Risk Factors Collaborators and others (2015).  Global, regional, and national comparative risk assessment of 79 behavioural, environmental and occupational, and metabolic risks or clusters of risks in 188 countries, 1990–2013: a systematic analysis for the global burden of disease study 2013.  Lancet 386(10010), 2287.
  * Haines et al. (1999) Haines, P. S., A. M. Siega-Riz, and B. M. Popkin (1999).  The diet quality index revised: a measurement instrument for populations.  Journal of the American Dietetic Association 99(6), 697–704.
  * Hashibe et al. (2007) Hashibe, M. et al. (2007).  Alcohol drinking in never users of tobacco, cigarette smoking in never drinkers, and the risk of head and neck cancer: pooled analysis in the international head and neck cancer epidemiology consortium.  Journal of the National Cancer Institute 99(10), 777–789.
  * Islam et al. (2019) Islam, M. A. et al. (2019).  Trans fatty acids and lipid profile: A serious risk factor to cardiovascular disease, cancer and diabetes.  Diabetes & Metabolic Syndrome: Clinical Research & Reviews 13(2), 1643–1647.
  * Johnson (1988) Johnson, R. A. R. A. (1988).  Applied multivariate statistical analysis (2nd ed. ed.).  Prentice-Hall series in statistics. Englewood Cliffs, N.J.: Prentice-Hall.
  * Kaiser (1958) Kaiser, H. F. (1958).  The varimax criterion for analytic rotation in factor analysis.  Psychometrika 23(3), 187–200.
  * Kromhout (2001) Kromhout, D. (2001).  Diet and cardiovascular diseases.  The journal of nutrition, health & aging 5(3), 144–149.
  * LaVange et al. (2010) LaVange, L. M. et al. (2010).  Sample design and cohort selection in the Hispanic Community Health Study/Study of Latinos.  Annals of epidemiology 20(8), 642–649.
  * Lichtenstein (2003) Lichtenstein, A. H. (2003).  Dietary fat and cardiovascular disease risk: quantity or quality?  Journal of Women’s Health 12(2), 109–114.
  * Liese et al. (2015) Liese, A. D. et al. (2015).  The dietary patterns methods project: synthesis of findings across cohorts and relevance to dietary guidance.  The Journal of nutrition 145(3), 393–402.
  * Liese et al. (2022) Liese, A. D. et al. (2022).  Variations in dietary patterns defined by the healthy eating index 2015 and associations with mortality: findings from the dietary patterns methods project.  The Journal of Nutrition 152(3), 796–804.
  * Lopes and West (2004) Lopes, H. F. and M. West (2004).  Bayesian model assessment in factor analysis.  Statistica Sinica 14(1), 41–67.
  * McLachlan (2008) McLachlan, G. J. (2008).  The EM algorithm and extensions (2nd ed. ed.).  Wiley series in probability and statistics. Hoboken, N.J.: Wiley-Interscience.
  * Meng and Rubin (1993) Meng, X.-L. and D. B. Rubin (1993).  Maximum likelihood estimation via the ecm algorithm: A general framework.  Biometrika 80(2), 267–278.
  * Morgenstern and Woodbury (1950) Morgenstern, O. and M. A. Woodbury (1950).  Stability of inverses of input-output matrices.  Econometrica 18, 190.
  * National Heart, Lung, and Blood Institute and others (2009) National Heart, Lung, and Blood Institute and others (2009).  Hispanic community health study/study of latinos (hchs/sol).  Bethesda, MD: NHLBI.
  * Robert and Escoufier (1976) Robert, P. and Y. Escoufier (1976).  A unifying tool for linear multivariate statistical methods: The rv- coefficient.  Applied Statistics 25(3), 257–265.
  * Schulze et al. (2003) Schulze, M. B., K. Hoffmann, A. Kroke, and H. Boeing (2003).  An approach to construct simplified measures of dietary patterns from exploratory factor analysis.  British Journal of Nutrition 89(3), 409–418.
  * Schwarz (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  The Annals of statistics 6(2), 461–464.
  * Sorlie et al. (2010) Sorlie, P. D. et al. (2010).  Design and implementation of the hispanic community health study/study of latinos.  Annals of epidemiology 20(8), 629–641.
  * Stearns et al. (2017) Stearns, J. C. et al. (2017).  Ethnic and diet-related differences in the healthy infant microbiome.  Genome medicine 9, 1–12.
  * Stephenson et al. (2020) Stephenson, B. J., A. H. Herring, and A. Olshan (2020).  Robust clustering with subpopulation-specific deviations.  Journal of the American Statistical Association 115(530), 521–537.
  * Trichopoulos et al. (1985) Trichopoulos, D. et al. (1985).  Diet and cancer of the stomach: a case-control study in greece.  International journal of cancer 36(3), 291–297.
  * Tucker et al. (1997) Tucker, L. A., G. T. Seljaas, and R. L. Hager (1997).  Body fat percentage of children varies according to their diet composition.  Journal of the American Dietetic Association 97(9), 981–986.
  * West (2003) West, M. (2003).  Bayesian factor regression models in the “large p, small n” paradigm, Volume 7.  Elsevier, North-Holland [u.a.].
  * Willett and MacMahon (1984) Willett, W. C. and B. MacMahon (1984).  Diet and cancer-an overview.  New England Journal of Medicine 310(11), 697–703.
  * Yoon et al. (2001) Yoon, P. W. et al. (2001).  The national birth defects prevention study.  Public health reports 116(Suppl 1), 32.
  * Zulyniak et al. (2017) Zulyniak, M. A. et al. (2017).  Does the impact of a plant-based diet during pregnancy on birth weight differ by ethnicity? a dietary pattern analysis from a prospective canadian birth cohort alliance.  BMJ open 7(11), e017753.