# Spatially Structured Regression for Non-conformable Spaces: Integrating
Pathology Imaging and Genomics Data in Cancer

Nathaniel Oshera,∗, Jian Kanga, Arvind Raob,c,a, Veerabhadran
Baladandayuthapania

aDepartment of Biostatistics, University of Michigan, Ann Arbor, MI 48105
bDepartment of Computational Medicine & Bioinformatics, University of
Michigan, Ann Arbor, MI
cDepartment of Radiation Oncology, University of Michigan, Ann Arbor, MI
∗Corresponding author. E-mail<EMAIL_ADDRESS>
(August 28, 2024)

###### Abstract

The spatial composition and cellular heterogeneity of the tumor
microenvironment plays a critical role in cancer development and progression.
High-definition pathology imaging of tumor biopsies provide a high-resolution
view of the spatial organization of different types of cells. This allows for
systematic assessment of intra- and inter-patient spatial cellular
interactions and heterogeneity by integrating accompanying patient-level
genomics data. However, joint modeling across tumor biopsies presents unique
challenges due to non-conformability (lack of a common spatial domain across
biopsies) as well as high-dimensionality. To address this problem, we propose
the Dual random effect and main effect selection model for Spatially
structured regression model (DreameSpase). DreameSpase employs a Bayesian
variable selection framework that facilitates the assessment of spatial
heterogeneity with respect to covariates both within (through fixed effects)
and between spaces (through spatial random effects) for non-conformable
spatial domains. We demonstrate the efficacy of DreameSpase via simulations
and integrative analyses of pathology imaging and gene expression data
obtained from $335$ melanoma biopsies. Our findings confirm several existing
relationships, e.g. neutrophil genes being associated with both inter- and
intra-patient spatial heterogeneity, as well as discovering novel
associations. We also provide freely available and computationally efficient
software for implementing DreameSpase.

Keywords: Bayesian variable selection, Cancer genomics, Cancer imaging, Data
integration, Random effect selection

##  1 Introduction

###  1.1 Scientific motivation and background

#### The tumor microenvironment

Tumors consist of many different components with complex interrelationships,
including the vascular network within the tumor, the extracellular matrix of
the tumor cells, and the immune cell composition of the tumor (Sun, 2016).
These structures are collectively referred to as the Tumor Microenvironment
(TME). The TME plays a critical role in the cancer development, progression
and the efficacy of certain treatments (Sadeghi Rad et al., 2021). Among its
various components, the immune composition of the TME has garnered much
interest among cancer researchers. This interest is due to increasing
recognition of cancer’s ability to avoid destruction by the immune system as
one of the organizing principles (“hallmarks”) of cancer evolution (Hanahan,
2022). Briefly, the immune composition of the TME is a fractally complex
ecosystem. The broad class of “immune cells” is divided into a multitude of
immune cell subtypes (Zoghi et al., 2023). Sometimes subtypes of the same
immune cell can even have divergent associations with prognosis (Li et al.,
2021). Adding to this complexity is the variability in the spatial composition
of the TME. While the prevalence of certain cell types within the TME can be
informative, more recent research has shown that the relative spatial
locations of different immune cell types (i.e. the tumor geography) offer
insights into tumor classification and prognostic assessment (Krishnan et al.,
2022).

#### Digital pathology

The complexity of TME has spurred research in the development and application
of systematic quantitative assessment and computational modeling, an area
often collectively referred to as Digital Pathology. Digital pathology
encompasses a broad array of aims and methods in the context of cancer
pathology imaging, including classification of cells in the TME, determination
of tumor grade, and the evaluation of spatial and morphological patterns in
the tumors (Heindl et al., 2015). Given the time and labor intensive nature of
assessing the spatial composition of the TME, the focus on evaluating spatial
and morphological patterns has become a central aim in digital pathology
research (Baxi et al., 2022). While there are many spatial characteristics of
tumors that are potentially of interest, such as the distributions of stromal
cells and blood vessels within the tumor, in this paper we focus on modeling
tumor-immune interaction. Our goal is to characterize the spatial
configurations of tumor cells and immune cells where the two cell types tend
to cluster as positively interacting, and configurations where they tend to
avoid one another as negatively interacting. This spatial tendency reflects
the infiltration of the tumor by immune cells.

#### Motivating dataset and key scientific questions

The motivating imaging data for this project comes from Melanoma biopsies
collected by The Cancer Imaging Archive (Clark et al., 2013). After some image
preprocessing, the cell-type (i.e. phenotypes) and locations of individual
cells in the biopsies are determined - see Figure 1 for examples of two such
biopsies (Osher et al., 2023); additional details provided in Section 5. In
addition to the high definition pathology images, bulk RNA-seq data were
collected for each patient at the biopsy level. This data naturally raises
three fundamental scientific questions. First, what is the appropriate method
of quantification of spatial tumor-immune interactions within and between
different biopsies (in other words, the extent to which immune cells interact
with tumor cells)? Second, how does the variability in these interactions
between biopsies relate to molecular composition of the tumor? Investigating
this may shed light on the biological mechanisms most responsible for tumor-
immune interactions. Third and finally, how does the variability in
interaction within biopsies (i.e. intra-tumoral heterogeneity) relate to
relate to molecular composition of the tumor? Fundamentally, answering these
three questions requires a synthesis of information from both imaging and
genomics data.

Figure 1: Plot of two biopsies from motivating Melanoma data set. Cell types
and locations for two biopsies in the data set. Black points represent tumor
cells, while red points represent immune cells. Enlarged portions are sub-
regions with varying degrees of tumor-immune interaction.

###  1.2 Statistical challenges

#### Data structure

As demonstrated by Figure 1, tumor-immune cell organization can vary
considerably both within and between biopsies, which can have drastically
different physical structures. Furthermore, it is of interest how these
cellular interactions vary both within and between biopsies with respect to a
set of covariates of interest (i.e. gene expression data). The resulting data
are complex and multifaceted, yet three key aspects are particularly
influential in determining the appropriate analysis approach. We refer to
these aspects as the resolution of the outcome, the resolution of the
covariates and the conformability of the spaces. The resolution of the outcome
refers to the “level” at which the outcome of interest is observed (i.e. at
the biopsy level or spatially across the biopsy), and the resolution of the
covariates is defined analogously. Conformability refers to whether or not the
spaces share a similar structure.

A great deal of work has been done in modeling of spatial data from
conformable spaces that can be “mapped” to a common space, along with high-
resolution outcomes and covariates (outcomes and covariates that are observed
at many spatial locations- see MacNab, 2022 and references there-in).
Additionally, a considerable amount of work has been done on data from
conformable spaces with either high-resolution outcomes and low-resolution
covariates or low-resolution outcomes and high-resolution covariates,
especially in the context of brain imaging data (Bowman, 2014 and references
there-in). However, comparatively little work has been done to model spatial
data from non-conformable spaces with high resolution outcomes and low
resolution covariates.

The structure of our motivating data set is considerably different than those
listed above. Each biopsy is partitioned into non-overlapping sub-regions,
which can then be modeled using spatial point processes. Measures of
interactions can then be summarized for each of the sub-regions, and modeled
with respect to biopsy level covariates such as gene expression data and/or
patient demographics and clinical outcomes. The result is data with high
resolution outcomes, low-resolution covariates, and non-conformable spaces.
While such structure poses no great obstacle to modeling heterogeneity with
respect to covariates between patients, modeling spatial heterogeneity within
patients with respect to covariates requires considerably more thought.

#### Existing methods

While there are existing methods that model interactions in the TME, they
largely focus on within biopsy modeling rather than joint modeling across
biopsies. Most notably in the statistical domain there has been work to
develop novel spatial models to assess interaction within the TME (Li et al.,
2019a ; Li et al., 2019b ). There have also been successful applications of
machine learning methods to this end, focusing on the application of deep
learning models and clustering (Abousamra et al., 2022; Saltz et al., 2018).
To our knowledge, a joint model to borrow strength across biopsies that would
answer the aforementioned scientific questions of interest has not been
developed. Another solution to this problem is to convert the non-conformable
spatial data to functional representations, which can then be analyzed in the
paradigm of functional data analysis (Vu et al., 2022; Chervoneva et al.,
2021; Yang et al., 2017). While this is a viable course of action for cross-
biopsy modeling, it does not take into account the inherent spatial structure
and correlation between adjacent sub-regions.

Figure 2: Overview of the DreameSpase model. (A) Point process modeling. Types
and locations are determined for each cell in high definition biopsy images
(top row to middle row). Biopsies are then divided into non-overlapping sub-
regions, on which a measure of spatial interaction is computed (middle row to
bottom row). (B) DreameSpase model fitting. Sub-regional outcomes are
vectorized and the DreameSpase model is fit on the resulting vectors. Gene
expression data, measured at the biopsy level, is used as the covariate set of
interest. (C) Fixed and Random Effect Selection. Selection of both fixed and
random effects are summarized for the model. fixed effects indicate
association between covariates and average level of infiltration in a given
biopsy, while random effects indicate association between covariates and
heterogeneity of infiltration in a given biopsy.

###  1.3 Model overview

#### Model components

In this paper, we present the Dual random effect and main effect selection
model for Spatially structured regression (DreameSpase). As conceptually
outlined in Figure 2, we first quantify the spatial interactions between
tumor-immune cells through point process modeling (Panel A). Subsequently,
DreameSpase model extends the spike and slab prior traditionally used for
fixed effects to the variance terms of spatial random effects. This allows for
the joint modeling of inter-biopsy heterogeneity with respect to covariates
via the fixed effects, and intra-biopsy heterogeneity with respect to
covariates via the spatial random effects (Panel B). The novelty of the model
lies in its ability to account for heterogeneity both within and between non-
conformable spaces with respect to arbitrary covariates of interest. The model
consists of two components of primary interest (Panel C):

  1. 1.

Fixed effect selection component. The model performs fixed effect selection
via spike and slab priors. In the context of cancer genomics, this component
accounts for the heterogeneity of the outcome of interest between biopsies
with respect to gene expression data.

  2. 2.

Spatial random effect selection component. In contrast to the fixed effect
selection component, the random effect selection component accounts for the
heterogeneity within biopsies with respect to gene expression data.

As we demonstrate via simulation studies, gains in selection accuracy through
joint modeling of the fixed and random effects are considerable as compared to
non-spatial and independent (two-step) approaches. We apply DreameSpase to our
motivating data set consisting of 335 images and associated gene expression
data for various immune classes. We find that the neutrophil class was most
consistently associated with within- and between-biopsy heterogeneity of
interaction between immune cells and tumor cells. The association between
several neutrophil class genes were consistent with the relationships
previously established in the literature, while others were suggestive of
potential novel relationships.

The remainder of the paper is organized as follows. Section 2 details the
formal model structure and formulation. Section 3 provides the details of the
sampling algorithm and inferential summaries of interest. Section 4 outlines
the simulation studies and assessments with comparative methods. Section 5
presents the results of analyzing whole slide Melanoma biopsy pathology
imaging and genomic data. Finally, we conclude with discussion, limitations,
and possible future directions for research in Section 6.

##  2 DreameSpase model

###  2.1 Modeling tumor-immune cell interactions

Figure 3: Illustration of interaction/co-localization. The leftmost panel
illustrates low-interaction, i.e. the cells of different types do not tend to
be near one another. The middle panel illustrates moderate interaction, i.e.
some cells of both types tend to be near one another, while others do not. The
third and final panel illustrates high interaction, i.e. the cells of both
types exhibit a strong tendency to be near cells of the other type.

#### The hierarchical Strauss model

In order to quantify the interaction between these two types of cells, we
employ a Bayesian Hierarchical Strauss Model (HSM) (Högmander and Särkkä,
1999). The HSM offers a parametric model that can capture both positive and
negative interaction between points of different types, i.e. the tendency of
the two types to either cluster towards one another or “avoid” one another. We
chose the HSM for two primary reasons. First, it can capture both positive and
negative interaction between different types of points, which the more
standard Strauss model cannot accomplish. Second, the HSM captures the
inherent biological nature of interaction between the tumor and the associated
immune cells. The HSM is called “hierarchical” because the modeling of the
second type of points (in our case, immune cells) is done conditionally on the
locations of the first type of points (in our case, tumor cells). While the
interaction between the tumor and immune cells is complex, as a first order
approximation it is natural to assume as a first order approximation that the
tumor precedes the immune response, and thus it makes intuitive sense to model
the immune response conditional on the locations of the tumor cells. The
density of the HSM model is as follows:

$f(\bm{x}_{1},\bm{x}_{2})\propto\exp\left\\{n_{1}\beta_{1}+n_{2}\beta_{2}+S_{R}(\bm{x}_{1},\bm{x}_{2})\theta\right\\},\vspace*{-10pt}$

where the parameters $\beta_{1}$ and $\beta_{2}$ capture the first order
intensity of points for types 1 and 2, respectively (tumor and immune cells),
with $n_{1}$ and $n_{2}$ representing the number of points for the
corresponding type. Of key importance in this model is the interaction
parameter, $\theta$, that captures the tendency of the different types of
points to exhibit positive or negative interaction. To gain some intuition as
to why this is the case, first note that $S_{R}(\bm{x}_{1},\bm{x}_{2})$ counts
the number of pairs of points within radius $R$ of one another where one is of
type 1 and the other is of type 2. Note that $R$ is not a parameter that is
fit, but is instead chosen prior to the model fitting based on knowledge of
the particular domain being modeled. Thus, when $\theta$ is negative it
“penalizing” the density put on configurations where points of the two types
tend to cluster. Conversely, such configurations being “rewarded” when
$\theta$ is positive. In essence, $\theta$ offers a highly interpretable
summary of the spatial interaction between the different types of points:

  * •

$\theta\in(-\infty,0)$ indicates negative interaction (as illustrated by the
leftmost panel of Figure 3), and lower values indicate more negative
interaction.

  * •

$\theta\in(0,\infty)$ indicates positive interaction (as illustrated by the
rightmost panel of Figure 3), and larger values indicate more positive
interaction.

  * •

$\theta=0$ indicates null interaction, i.e. the points are “indifferent” to
points of the other type, as illustrated by the middle panel of Figure 3.

#### Biopsy partitioning

Whole slide high-resolution pathology images of biopsies typically contain
between $10^{4}-10^{5}$ cells. For example, in our motivating data set, the
median number of cells is $\sim 50,000$, with the maximum across all biopsies
$\sim 190,000$. This volume means that simply computing a spatial summary on
the biopsy level is challenging. Due to the sheer size of the biopsy any given
spatial feature will tend to vary across the biopsy in ways that may be
important. We are primarily interested in modeling interaction (the tendency
of different cells to be near each other), which is quantified on a fairly
local scale ($\sim 30\mu m$) which is dictated by radius of influence of a
single cell. In order to capture the heterogeneity of this interaction, we
partition the biopsy into sub-regions defined by the locations of the
individual tumor cells. In addition, this partitioning has the added benefit
of improving the computational requirements of the HSM.

More specifically, let $N$ be the number of biopsies, respectively denoted
$B_{1},\dots,B_{N}$, each $B_{i}\subset\mathbb{R}^{2}$ is further partitioned
into $n_{i}$ disjoint sub-regions $b_{i(1)},\dots,b_{i(n_{i})}$ with
$b_{i(j)}\subseteq B_{i}$, $b_{i(j)}\cap b_{i(k)}=\varnothing$ for $k\neq j$,
and $B_{i}=\bigcup\limits_{j=1}^{n_{i}}b_{i(j)}$. On sub-region $b_{i(j)}$ we
observe a marked point process model $\bm{P}_{i(j)}$ consisting of two types
of points, tumor cell and immune cell.

The HSM can be fit on each of the non-overlapping sub-regions to determine the
local level of tumor-immune interaction, as shown in panel A of Figure 2. For
more details on the HSM and the $\theta$ parameter, see Sections S1.2 and S1.3
of the supplementary materials. The model can be fit using standard Bayesian
posterior sampling methods on the so-called pseudolikelihood function used in
the estimation of these models; see Section S1.2 of the supplementary
materials for model fitting details. While there are many options in the
spatial statistical literature for measuring and summarizing interaction
between different types of points, one of the primary advantages of $\theta$
is that it tends to be normally distributed in its posterior distribution at
the level of sub-regions (Figure S1 in the supplementary materials).

At the end of this modeling step, we obtain a partition of each biopsy into
non-overlapping sub-regions, and a measure of tumor-immune interaction on each
of the sub-regions. Let $N$ be the number of biopsies, denoted
$B_{1},\dots,B_{N}$, where biopsy $B_{i}$ has $n_{i}$ observed cells
$c_{i1},\dots,c_{i{n_{i}}}$. For each $c_{ij}$, two pieces of information are
observed: its location on the slide, $\bm{x}_{ij}=[x_{ij1},x_{ij2}]$, and its
type, $t_{ij}$, which is either “tumor” or “immune.” Next, we specify a joint
model for the tumor-immune interaction across biopsies conditional on biopsy
level covariates.

###  2.2 Spatially structured regression model

From the previous step the resulting data structure is as follows: for each
biopsy $B_{1},\dots,B_{N}$, suppose we observe on each biopsy $B_{i}$ a
continuous outcome of interest at locations $\bm{S}_{i}$, i.e.
$\bm{Y}_{i}(\bm{S}_{i})=[y_{i}(s_{i(1)}),\dots,y_{i}(s_{i}(n_{i}))]^{T}$. Note
that in our application, $y_{i}(s_{i(1)})$ is the $\theta$ parameter estimated
on that sub-region as described in Section 2.1. Additionally, suppose that for
each biopsy, $p$ biopsy-level covariates (in our application, gene expression
data) are observed, denoted by $\bm{X}_{i}=[x_{i1},\dots,x_{ip}]^{T}$. The
primary motivation for our modeling is to assess whether or not the biopsy
level covariates are associated at the population level with increased
heterogeneity within or between biopsies. To this end, we define the
DreameSpase model for biopsy $i\in\\{1,\ldots,N\\}$ as:

$\bm{Y}_{i}(\bm{S}_{i})=(\overbrace{\bm{\alpha}^{T}}^{\text{Fixed
effects}}\bm{X}_{i})+\overbrace{\bm{\eta}_{i}(\bm{S}_{i})}^{\text{Covariate
specific random
effects}}\bm{X}_{i}+\overbrace{\bm{\delta}_{i}(\bm{S}_{i})}^{\text{Spatial
random effect}}+\bm{\epsilon}_{i}.$ (1)

We discuss each component of the model in detail below.

#### Fixed effects

In (1), $\bm{\alpha}^{T}\bm{X}_{i}$ is the biopsy level mean of sub-region
level outcomes. Thus, $\bm{\alpha}$ captures the variance explained by the
covariates between different biopsies and, conditional upon the covariates,
can be interpreted in expectation at a population level much like the
coefficients of a standard multiple regression. The object of this component
of the model is to identify covariates which are systematically associated
with an increase the average level of interaction in a given biopsy.

#### Covariate specific random effects

The term $\bm{\eta}_{i}(\bm{S}_{i})$ in (1) is an $n_{i}\times p$ matrix of
covariate specific spatial random effects for each biopsy. Each column of
$\bm{\eta}_{i}(\bm{S}_{i})$ represents a spatial random effect for a specific
covariate. More concretely, denoting the $j$th column of $\bm{\eta}_{i}$ by
$\bm{\eta}_{i(j)}(\bm{S}_{i})$, we assume that $\bm{\eta}_{i(j)}(\bm{S}_{i})$
follows a conditionally autoregressive (CAR) process. The CAR process and
other models derived from it are commonly used to analyze areal data in order
to share information across adjacent sub-regions (Orozco-Acosta et al., 2022).
The form of the specific CAR process is given by:

$\bm{\eta}_{i(j)}(\bm{S}_{i})\thicksim
MVN(\bm{0},\psi_{j}^{2}\left[D_{w}(\bm{W}_{i})-\phi\bm{W}_{i}\right]^{-1}).$
(2)

While other options for modeling areal data exist, the CAR prior is appealing
since it makes minimal assumptions about the structure of the spatial
correlation beyond the dependency of regions on their neighbors as well as the
fairly small number of parameters. This is particularly appealing given the
non-conformability of the underlying tumor biopsies. Under this specification,
$\bm{W}_{i}$ is an $n_{i}\times n_{i}$ symmetric matrix, where the entry in
column $j$ of row $i$ (and column $i$ of row $j$) is 1 if sub-regions $i$ and
$j$ are adjacent, and 0 otherwise; $D_{w}(\bm{W}_{i})$ is a diagonal matrix
where each diagonal entry is the sum of the corresponding row of $\bm{W}_{i}$;
$\psi_{j}^{2}$ captures the variability of the value in a sub-region
conditional upon the value of its neighbors; and $\phi$ captures the spatial
correlation between neighbors. Note that while $\phi$ takes on values in
$(-1,1)$ and determines the correlation between adjacent sub-regions under
this model specification, it is not itself equal to the correlation between
adjacent sub-regions. The distribution of this spatial random effect is thus
determined by the values of $\psi_{j}^{2}$ and $\phi$.

The key parameter when considering the covariate specific spatial random
effects is $\psi_{j}^{2}$. This parameter affects the outcome via the
multiplication between the matrix of random effects and the covariate vector.
Consider, as an extreme but illustrative example, the case where
$\psi_{j}^{2}=0$. When this is the case, we should expect no systematic change
heterogeneity within a biopsy as the $j$th covariate changes, as each
$\bm{\eta}_{i(j)}(\bm{S}_{i})$ has a degenerate distribution at $\bm{0}$.
However, consider $\psi_{j}^{2}>0$. As the magnitude of $x_{ij}$ increases, we
should in turn expect to see a systematic increase in the heterogeneity of
biopsy $i$. Note that this may not translate to an increase in the average
level of heterogeneity in biopsy $i$\- whether or not such an increase occurs
depends on the value of $\alpha_{j}$. These two examples, while extreme,
embody the logic underlying our decision to use spike and slab priors to
regularize the fittings of these parameters; details deferred to Section 2.3
for now.

#### Global CAR process

In (1), $\bm{\delta}_{i}(\bm{S}_{i})$ is the “global” spatial random effect
that captures the overall spatial correlation across all biopsies. We assume
that that $\bm{\delta}_{i}(\bm{S}_{i})$ follows an analogous CAR process,
defined as:

$\bm{\delta}_{i}(\bm{S}_{i})\thicksim
MVN(\bm{0},\tau^{2}\left[D_{w}(\bm{W}_{i})-\rho\bm{W}_{i}\right]^{-1}),$

where each term is defined analogously to its counterpart in (2). Like the
covariate specific CAR process, the distribution of this spatial random effect
is thus determined by the values of $\tau^{2}$ and $\rho$. However, unlike the
covariate specific CAR process, the $\rho$ parameter is estimated. While the
model is well-defined for $\rho\in(-1,1)$, values less than zero indicate
anti-correlation between neighbors. Because this is not a priori biologically
plausible (there is no reason to suspect that adjacent tumor sub-regions
should differ systematically), we assume that $\rho\in[0,1)$, i.e. that sub-
regions are on average either positively correlated or independent of their
neighbors.

#### Pure error

Finally, the $\bm{\epsilon}_{i}\thicksim MVN(\bm{0},\nu^{2}\bm{I})$ is the
pure error of the model (often referred to as the “nugget”), i.e. the error
that is not accounted for by the fixed effects, the global CAR process, or the
covariate specific CAR processes.

It should be noted that the DreameSpase model is agnostic to how to the values
on the sub-regions are generated; nothing about the downstream model depends
on the usage of point processes to generate sub-regional outcome measurements.
The only assumption is that the outcomes observed at the level of each sub-
region are spatially correlated and (approximately) normally distributed.

###  2.3 Dual main and random effect spatial selection

#### Motivation

A challenge of working with genomic data is the dimensionality. Even when
utilizing a targeted subset of available genes, it is neither likely nor
expected that all genes included in the model are meaningfully related to the
outcome. This necessitates some degree of selection or regularization in order
to properly determine which genes are significantly associated with the
outcome of interest. Performing selection yields the dual benefits of
producing a more parsimonious model while also improving the accuracy of
parameter estimates. To accomplish this in the context of our model, we
utilize spike and slab priors for both the fixed and random effects. We focus
discussion on these priors due to their importance in the model.

#### Spatial fixed effects

While there are several ways to parameterize a spike and slab prior, our
method represents a standard application of stochastic search variable
selection (George and McCulloch, (1993)). The prior is specified for each
entry of the fixed effect vector $\alpha_{j}$ as follows:

$\alpha_{j}\thicksim\gamma_{j}N(0,\sigma^{2}_{slab})+(1-\gamma_{j})N(0,\sigma^{2}_{spike}),$
$\gamma_{j}\thicksim Bernoulli(P_{\gamma}),$

where $P_{\gamma}$ is the prior probability of inclusion for the fixed effects
and is assumed to come from a common (non-informative) $Beta(1,1)$ prior
distribution. Note that under our specification, as is often the case,
$\sigma^{2}_{spike}$ and $\sigma^{2}_{slab}$ are not estimated, but rather
selected a priori; in our application, we set $\sigma^{2}_{spike}=0.03$ and
$\sigma^{2}_{slab}=100$. This prior construction allows not only the
estimation of each $\alpha_{j}$ via the posterior mean, but also the
probability of inclusion in the model via the posterior mean of the
$\gamma_{j}$ parameter. Because $\gamma_{j}$ is an indicator, the posterior
expectation corresponds to the probability that $\gamma_{j}=1$, which is the
probability that $\alpha_{j}$ is included in the model (i.e. is not penalized
by the prior).

#### Spatial random effects

Unlike standard fixed effect coefficients, variance terms are constrained by
definition to lie in the range $[0,\infty)$. The challenge in adapting the
spike and slab prior to the variance terms of random effects stems from the
difference in domains between the variance terms
$\psi_{1}^{2},\dots,\psi_{p}^{2}$ and the fixed effect terms
$\alpha_{1},\dots,\alpha_{p}$. In theory each fixed effect term can take on
any value in $\mathbb{R}$, even if excessively large or small values are
implausible in practice. This is not the case when dealing with variance
parameters. By definition, such parameters must lie in $[0,\infty)$. This
rules out the usage of normal priors centered at zero, which put half their
mass in the range $(-\infty,0]$. Because of this asymmetry, various models
have been constructed by utilizing specialized Cholesky decompositions of the
covariance matrix of the random effect terms (Joyner et al., 2020; Ibrahim et
al., 2011). By placing penalties (in the form of priors or regularization
components) on terms in the decompositions of the covariance matrix rather
than directly on the terms of the covariance matrix, these models avoid
dealing with both the constraints of covariance matrices as well as the domain
restrictions of variance terms. As a result, such penalties can be constructed
in the same manner as they might for more standard fixed effect terms. Others,
such as Scheipl et al., 2012, employ a spike and slab formulation on variance
terms in the context of structural additive regression (STAR) models. However,
the form and purpose of this model differs considerably from our own, as the
purpose of STAR models is to provide a flexible and generalizable method for
modeling the mean of some outcome vector.

We define the “half-spike and half-slab” priors for the variance terms using
half-normal distributions as follows:

$\psi^{2}_{j}\thicksim
d_{j}N^{+}(0,\xi^{2}_{slab})+(1-d_{j})N^{+}(0,\xi^{2}_{spike}),$
$d_{j}\thicksim Bernoulli(P_{d}),$

where $P_{d}\thicksim Beta(1,1)$ is the prior probability of inclusion for all
random effects. Like those of the fixed effects, the half spike is
parameterized by a small variance, $\xi^{2}_{spike}$, while the slab is
parameterized by a larger variance, $\xi^{2}_{slab}$. Half-normal
distributions were used in place of more standard inverse gamma distributions
partially due to the relatively larger amount of mass placed on extreme values
by the half-normal and to achieve greater congruence between the fixed effect
selection component of the model and random effect selection component of the
model. This congruence makes it more straightforward to tune the selection of
both the fixed and random effects in the model. The $d_{j}$ indicator variable
also serves a similar function to the $\gamma_{j}$ in the fixed effect
formulation.

Table 1: Full specification of DreameSpase model and priors. Model specification is given for $N$ biopsies with $p$ covariates of interest. Across all terms, $i$ indexes the biopsy and $j$ indexes the covariate where relevant. $\bm{\eta}_{i(j)}$ is the $j$th column of the $n_{i}\times p$ matrix $\bm{\eta}_{i}$, the matrix of covariate specific random effects for biopsy $i$. $\bm{Y}_{i}\mid\bm{\mu}_{i},\nu^{2}$ | $\thicksim$ |  $MVN(\bm{\mu}_{i},\nu^{2}\bm{I})$, | $i=1,\dots,N$
---|---|---|---
$\bm{\mu}_{i}$ | $=$ | $\left(\bm{\alpha}^{T}\bm{X}_{i}\right)\cdot\bm{1}+\bm{\delta}_{i}(\bm{S}_{i})+\bm{\eta}_{i}(\bm{S}_{i})\bm{X}_{i}$ | $i=1,\dots,N$
$\alpha_{j}\mid\gamma_{j},\sigma^{2}_{spike},\sigma^{2}_{slab}$ | $\thicksim$ |  $\gamma_{j}N(0,\sigma^{2}_{slab})+(1-\gamma_{j})N(0,\sigma^{2}_{spike})$, | $j=1,\dots,p$
$\bm{\delta}_{i}(\bm{S}_{i})\mid\tau^{2},\rho,\bm{W}_{i}$ | $\thicksim$ |  $MVN(\bm{0},\tau^{2}\left[D_{w}(\bm{W}_{i})-\rho\bm{W}_{i}\right]^{-1})$, | $i=1,\dots,N$
$\bm{\eta}_{i(j)}(\bm{S}_{i})\mid\psi_{j},\phi,\bm{W}_{i}$ | $\thicksim$ | $MVN(\bm{0},\psi_{j}^{2}\left[D_{w}(\bm{W}_{i})-\phi\bm{W}_{i}\right]^{-1})$ |  $i=1,\dots,N$; $j=1,\dots,p$
$\psi_{j}\mid d_{j},\xi^{2}_{spike},\xi^{2}_{slab}$ | $\thicksim$ | $d_{j}N^{+}(0,\xi^{2}_{slab})+(1-d_{j})N^{+}(0,\xi^{2}_{spike})$ | $j=1,\dots,p$
$\nu^{2}$, $\tau^{2}$ | $\thicksim$ | $Inverse\leavevmode\nobreak\ Gamma(0.001,0.001)$ |
$\rho$ | $\thicksim$ | $Discrete(v_{1},\dots,v_{K},p_{1},\dots,p_{K})$ |  $v_{\ell}\in[0,1)$; $p_{\ell}\in[0,1]$; $\sum_{\ell=1}^{K}p_{\ell}=1$
$\gamma_{j}\mid P_{\gamma}$ | $\thicksim$ |  $Bernoulli(P_{\gamma})$, | $j=1,\dots,p$
$d_{j}\mid P_{d}$ | $\thicksim$ |  $Bernoulli(P_{d})$, | $j=1,\dots,p$
$P_{\gamma},P_{d}$ | $\thicksim$ | $Beta(1,1)$ |

#### Additional priors

The remaining priors in the model are generally chosen to be standard
conjugate priors, with the exception of the $\rho$ parameter. $\rho$ is given
a discrete distribution $P(\rho=v_{i})=p_{v_{i}}$ for $i\in\\{1,\dots,K\\}$;
see Section 3.1 for details. As a brief point of discussion, we have chosen a
non-informative $Beta(1,1)$ prior for the prior selection probabilities for
the fixed and random effects, this prior could instead be assigned uniquely to
each selection probability and parameterized based on prior clinical
information on the importance of each covariate in the style of Ni et al.,
2019. For full prior specifications, see Table 1.

##  3 Posterior inference

###  3.1 Key full conditional distributions

We perform posterior inference via Markov Chain Monte Carlo on the full
parameter set outlined above: fixed effects $\bm{\alpha}$, corresponding
selection indicators $\bm{\gamma}=[\gamma_{1},\dots,\gamma_{p}]^{T}$, and
overall probability of selection $P_{\gamma}$; random effect variance terms
$\bm{\psi}=[\psi_{1}^{2},\dots,\psi_{p}^{2}]$, corresponding selection
indicators $\bm{d}=[d_{1},\dots,d_{p}]^{T}$, and overall probability of
selection $P_{d}$; the global CAR process parameters $\tau^{2}$ and $\rho$;
and the pure error variance $\nu^{2}$. In this section we outline the
derivations of selected full conditional distributions necessary for the
implementation of the Gibbs sampling algorithm. While many of the full
conditional distributions utilize conjugate or otherwise standard priors,
there are several prior specifications that result in full conditional
distributions unique to our model. See Table 1 for the full model
specification from which all full conditional distributions can be derived,
and Section S3.3 for detailed derivations of all model parameters. Due to
identifiability constraints, $\phi$ is fixed a priori at a relatively small
value (0.3).

#### Spatial random effects

The full conditional of $\psi_{j}^{2}$ is induced via the prior on $\psi_{j}$,
and requires somewhat more careful derivation. Specifying a half-normal prior
on $\psi_{j}$ causes the full conditional distribution of $\psi_{j}^{2}$ to be
that of a generalized inverse Gaussian distribution. For a detailed derivation
of this fact, see Section S3.3.2 of the supplementary materials. Thus, while
the prior is specified in terms of the standard deviation $\psi_{j}$, in
practice the sampling is performed on the variance $\psi_{j}^{2}$. However,
due to the simple monotonic relationship between the two quantities, this does
not pose any practical inconvenience.

#### Global CAR process correlations

In order to improve the computational efficiency of sampling, a grid prior was
placed on the correlation of the global CAR process, $\rho$. This prior is
specified by a discrete set of values and corresponding probabilities,
$v_{1},\dots,v_{K}$, and $p_{1},\dots,p_{K}$ (full conditional distributions
provided in Supplementary Section S3.3.3). In practice, sampling from this
distribution is accomplished using the Gumbel Max trick (Huijben et al.,
2023); see Supplementary Section S3.1 details and derivation.

###  3.2 Inferential summaries of interest

There are two primary inferential summaries of interest: the sub-region level
summaries of interaction which are summarized within and across biopsies, and
the selection probabilities for the fixed effects and random effects. While we
are primarily interested in the former only as the object of modeling in this
paper, these summaries can be explored in many ways in their own right. As for
the latter, the usage of spike and slab priors allows for the computation of
posterior probabilities of selection for both the fixed and random effects
$\bm{d}$ and $\bm{\Gamma}$, illustrated in panel C of Figure 2. For each
covariate, any combination of the corresponding fixed and random effect can be
selected for inclusion in the model as explained in Section 2.3. In the
context of our application:

  * •

$P(\gamma_{j}=1)$ is the probability that the $j$th fixed effect is included
in the model, indicating that the corresponding gene is associated with
heterogeneity between biopsies

  * •

$P(d_{j}=1)$ is the probability that the $j$th random effect is included in
the model, indicating that the corresponding gene is associated with
heterogeneity within biopsies

If both the fixed and random effects are selected, then a given gene is
associated with both inter- and intra-biopsy heterogeneity, and thus may be of
particular interest.

##  4 Simulation studies

#### Simulation overview

In order to assess the overall performance of our DreameSpase model we
performed a series of simulation studies. In terms of design parameters, the
number of biopsies $N$ was fixed at 200 (on the same order as a real data
example), while the number of covariates $p$ was varied between 50, 75, and 90
to achieve ratios of $2p/N$ of 0.5, 0.75, and 0.9. We refer to the ratio
$2p/N$ as the relative dimensionality of a setting. All covariates were
generated from a standard normal distribution. The number of true fixed and
random effects for each of the three settings was fixed at 3, 6, and 9,
respectively. Each biopsy was simulated as a $5\times 5$ lattice across all
settings. Fixed and random effects were simulated with equal frequency to be
small, medium, or large. Small, medium, and large fixed effects were simulated
uniformly in the respective ranges $(0.13,0.23)$, $(0.23,0.33)$, and
$(0.33,0.43)$, while the corresponding ranges for random effect variance were
$(0.1,0.2)$, $(0.2,0.3)$, and $(0.3,0.4)$. 50 data sets were simulated for
each data setting.

#### Comparative methods

In addition to the DreameSpase model, two additional models were fit to serve
as performance benchmarks. The first model is a special case of the
DreameSpase model that uses non-spatial random effects. This is equivalent to
setting $\phi=0$ and $\bm{W}_{i}=\mathrm{diag}(1,\ldots,1)$ in Equation 2.
This ignores the spatial structure of the biopsy and instead models the
$\bm{\eta}_{i(j)}$ as simply clustered data. As in the true model effects were
selected if their estimated posterior probability of inclusion was greater
than $0.5$ in accordance with the median-probability selection model (Barbieri
and Berger, 2004). We refer to this model as the “Non-Spatial DreameSpase”
model, abbreviated to NSDS.

The second method was a penalization-based method that a researcher might use
to answer the questions posed in this paper without developing a novel method,
which we refer to as the “analyst model.” There were two stages of analysis
for this model: one for the fixed effect component, and one for the random
effect selection component. For the fixed effect selection component, a lasso
regression was fit using the sub-region level outcomes as the outcome and the
biopsy level covariates as the covariates. A penalization parameter was
selected via three-fold cross validation, and fixed effects were considered to
be selected if their parameter estimates were non-zero. To perform random
effect selection, a lasso regression was fitted treating the standard
deviation of each biopsy as the outcome, and the absolute value of the
covariates as the regressors. In this case as well the penalization parameter
was chosen via three-fold cross validation, and non-zero covariates were
considered to be selected.

#### Evaluation criteria

The primary goal was to assess selection accuracy of the different models and
methods. To do so, the true positive rate and false positive rate of selection
were computed for the fixed effects as well as the random effects. The true
positive rate is computed by the size of the effects, while the false positive
rate is presented across all effect sizes (since null effects do not have
sizes).

#### Signal to noise ratio

In order to properly calibrate the simulation studies to the real data and
compare across differing settings, it is useful to define a signal to noise
ratio. Note that unlike many models, we are interested in accurately
estimating not only the fixed effect terms, but random effect terms as well.
Because the two components of the model are independent of one another, we
define separate signal to noise ratios for both components of the model. The
fixed and random effect SNRs are defined as

$SNR_{fixed}=\frac{\text{var}(\bm{X}^{T}\bm{\alpha})}{\text{var}(\bm{Y})},\qquad
SNR_{rand}=\frac{\sum_{j=1}^{n}\text{var}(\bm{\eta}_{j,\cdot}\bm{X})}{n\text{var}(\bm{Y})}.$

See Section S2.1 for details and derivations for both quantities. The chosen
simulation settings yielded simulations with fixed effect SNRs in the range
$[0.37,0.42]$, and random effect SNRs in the range $[0.38,0.39]$.

#### Results

Table 2 shows the detailed results of the true positive rates across the
effect sizes and relative dimensionality and Table 3 shows for results of the
false positive rates.

Across all settings, the false positive rate of the DreameSpase model never
exceeded $0.02$. For medium and large effect sizes, the true positive rates
for the fixed effects are uniformly $1.0$, i.e. perfect selection. For the
random effects, the true positive rate lies between 0.8 and 1.0, with the
median for these effect sizes being $0.98$. The model struggles substantially
more with smaller effect sizes, with the lowest true positive rate obtained
being 0.44 when $2p/N=0.9$. Still, this is only marginally worse than the
analyst model, which achieved the best true positive rate for this setting at
$0.53$. Thus, while the model may not be sufficiently sensitive to detect all
effects of this magnitude, the highly controlled false positive rate means
that one can be confident that the effects we do find are true effects.
Moreover, while DreameSpase is sometimes outperformed in selection of the
smallest magnitude effects by the other methods, the differences in
performance are generally quite small, while the deficiencies of the other
methods across the different effect types and sizes are generally quite
substantial.

The analyst model was able to select random effects with reasonable accuracy,
and in fact outperforms DreameSpase as well as the NSDS model in the
$2p/N=0.9$ setting in terms of true positive rate, consistently achieving the
best selection across several settings. The lowest true positive rate for
random effect selection achieved is 0.56, (which is still the best across all
three models at this setting), while the highest is 1.0. The false positive
rate is also well controlled for random effects, again never exceeding 0.02.

However, the failure to jointly model fixed effects and random effects
resulted in a severely diminished true positive rate and inflated false
positive rate for the fixed effects not seen when using the DreameSpase model.
Because of the failure to jointly model the different effect types, the true
positive rate for fixed effects lay between 0.4 and 0.56. Additionally, the
false positive rate lay between 0.21 and 0.24, while DreameSpase and NSDS
achieved false positive rates of less than 0.01. The only exception to this
trend is the true positive rate among small effects in the $2p/N=0.9$ setting,
where the analyst model achieves a true positive rate of $0.53$ compared to
0.44 for both the DreameSpase as well as the NSDS model.

Table 2: True positive rate (TPR) across settings, models, and effect sizes. Relative dimensionality is defined as $2p/N$, where $N$ is the number of biopsies and $p$ is the number of covariates. The number of covariates is doubled in the ratio to account for the estimation of both fixed effects and covariate specific random effects, thus doubling the effective dimensionality. Effect sizes refer to the magnitude of the fixed and random effects. For exact ranges, see Section 4. Bolded values indicate the best performance for a given setting and effect size. In the case of a tie, all values are bolded. TPR |  | Fixed Effects | Random Effects
---|---|---|---
$2p/N$ | Effect Size | “Analyst Model” | NSDSE | DreameSpase | “Analyst Model” | NSDSE | DreameSpase
| Small | 0.56 | 1.00 | 1.00 | 0.94 | 0.98 | 0.98
0.5 | Medium | 0.46 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
| Large | 0.4 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
| Small | 0.48 | 0.63 | 0.61 | 0.76 | 0.26 | 0.76
0.75 | Medium | 0.54 | 1.00 | 1.00 | 0.99 | 0.75 | 0.96
| Large | 0.51 | 1.00 | 1.00 | 1.00 | 0.93 | 1.00
| Small | 0.53 | 0.44 | 0.44 | 0.56 | 0.05 | 0.53
0.9 | Medium | 0.51 | 1.00 | 1.00 | 0.93 | 0.12 | 0.8
| Large | 0.49 | 1.00 | 1.00 | 0.99 | 0.2 | 0.93

For fixed effects, the performance of the NSDS model is comparable to that of
DreameSpase across all settings in terms of both TPR and FPR. For random
effects, however, the true positive rate of the NSDS model is similar to that
of DreameSpase in the $2p/N=0.5$ setting, but deteriorates quickly across the
larger settings, at times reaching as low as 0.05. Notably, in the $2p/N=0.9$
setting, the maximum true positive rate for random effects achieved by the
NSDS model across all effect sizes is 0.2, compared to 0.93 for the
DreameSpase model and 0.99 for the analyst model. While it achieves a false
positive rate of virtually 0 in this setting, the DreameSpase and analyst
models both achieve a modest false positive rate while also achieving a
substantially higher true positive rate across all effect sizes. This
indicates the importance of accounting for the spatial structure of the
biopsies for the random effects.

#### Additional simulation results

We also computed threshold free metrics, in particular Area under the curve
(AUC) for all the methods and the results are shown in Table S1 of the
supplementary materials. The results are broadly consistent with the ones
presented above; the analyst model tends to perform well in random effect
selection (normalized AUC ranging from 0.92 to 1), but struggles to a
significant degree in main effect selection (normalized AUC ranging from 0.47
to 0.51).

Table 3: False positive rate (FPR) across settings and models. Relative dimensionality is defined as $2p/N$, where $N$ is the number of biopsies and $p$ is the number of covariates. The number of covariates is doubled in the ratio to account for the estimation of both fixed effects and covariate specific random effects, thus doubling the effective dimensionality. Note that false positive rates are determined solely by the selection results for the simulated null covariates, and thus the FPR does not vary when broken down by effect size. Bolded values indicate the best performance for a given setting and effect size. In the case of a tie, all values are bolded. FPR | Fixed Effects | Random Effects
---|---|---
$2p/N$ | “Analyst Model” | NSDS | DreameSpase | “Analyst Model” | NSDS | DreameSpase
0.5 | 0.24 | < 0.01 | < 0.01 | < 0.01 | < 0.01 | 0.02
0.75 | 0.22 | < 0.01 | < 0.01 | 0.01 | < 0.01 | 0.01
0.9 | 0.21 | < 0.01 | < 0.01 | 0.02 | < 0.01 | 0.01

##  5 Tumor microenvironment in melanoma

###  5.1 Background

#### Melanoma and immunotherapy

Melanoma is a particularly interesting type of cancer given its high
immunogenicity caused by its high mutational load (Marzagalli et al., 2019).
This characteristic has implications not only for patient prognosis, but also
susceptibility to treatment via immunotherapies (Kang et al., 2020; Simiczyjew
et al., 2020). While remarkably effective in certain patients, immunotherapy
is by no means a silver bullet in the treatment of melanoma. Ipilimumab, the
first CTLA-4 inhibitor approved by the FDA for treatment of melanoma patients,
has demonstrated remarkable efficacy in the treatment of even advanced stage
melanoma, both alone and in concert with other treatments. However, the
objective response rate across these different settings is far from 100%
across all patients (Chesney et al., 2018). Further complicating this picture
is the fact that the composition of TME has also been shown to be associated
with response to immunotherapies (Falcone et al., 2020). The biological
understanding of the interplay between immunotherapy and the composition of
the TME has been bolstered recently by the proliferation of single cell
imaging technologies (Xiao et al., 2022). Still, more work is required to
synthesize our understandings of immunological mechanisms in melanoma and its
relationship to the spatial composition of the TME. To this end, we
investigate the relationship between spatial composition of the TME (using
pathology imaging data) and immunological mechanisms (using genomic data),
both at the level of tumor-immune interactions as well as the heterogeneity of
the interactions within and across tumors.

#### Description of pathology imaging and genomic data

Our motivating imaging dataset comes from the The Cancer Imaging Archive
(Clark et al., 2013), a companion consortium to the The Cancer Genome Atlas
project (TCGA). TCGA is an initiative that makes available high quality
genomics and imaging data on samples from a diverse range of cancer types
(Weinstein et al., 2013). This multi-platform data includes, most notably for
our purposes, high definition images of the biopsies and gene expression data
from patient samples. Our pathology imaging data set consisted of whole slide
imaging data for 335 biopsies from the Skin Cutaneous Melanoma (SKCM), with
the number of cells per biopsy ranging from 814 to 187,521. The cells of each
biopsy were processed and classified by a machine learning model as either
tumor or immune according to the morphology of each cell. Biopsies were
subsequently partitioned into non-overlapping sub-regions, and an estimate of
interaction was computed for each sub-region. For more details on the cell
classification algorithm and tumor partitioning algorithm, see Osher et al.,
(2023).

In order to determine the association between immunological mechanisms and
spatial association of tumor cells and immune cells, we used gene expression
data from genes associated with different groups of immune cells. These immune
groups include “marker” genes linked to B cells, macrophages, monocytes,
neutrophils, natural killer cells, plasma cells, and T cells. These immune
cell groups regulate many important oncogenic processes such as angiogenesis,
metastasis, and the response to drugs or immunotherapy (Nirmal et al., 2018).

Using the sub-regional interaction parameter estimates described in Section
2.1 as the outcome, the DreameSpase model was fit using gene expression data
for the previously described sets of genes downloaded from the UCSC Xena
browser (Goldman et al., 2020). In total, there were 11,202 sub-regions on
which interaction was estimated across 335 biopsies, with the median number of
sub-regions per biopsy being 29 (IQR 26). After pre-processing, there were 122
genes included in our final analysis and additional pre-processing details are
available in Section S3.2 of the supplementary materials; see Algorithm S1 for
the exact algorithm used.

The DreameSpase model was run for 300,000 total iterations, 150,000 of which
were treated as burn-in. The resulting samples were thinned by 10, resulting
in 15,000 posterior samples. Global convergence for the posterior samples was
assessed via the Geweke diagnostic on the likelihood ($p>0.05$); see Figure S2
of the supplementary materials for trace plots of log-likelihood and selected
parameters. The analysis took approximately 8 hours on an M1 Macbook Pro with
8 gigabytes of memory.

###  5.2 Results

#### Overview

Figure 4 illustrates the spatial heterogeneity of the tumor-immune interaction
both within biopsies (captured by the standard deviation of sub-regional
$\theta$ estimates) and between them (captured by the mean of the $\theta$
estimates). Biopsy level means ranged from -0.09 to 0.49, while biopsy level
standard deviations ranged from 0.02 to 0.32. The absolute value of the
coefficients of variation also ranged considerably, the minimum being 0.08 and
the maximum being 437. This illustrates the degree to which the average level
of interaction within a biopsy can vary separately from the level of
heterogeneity within a biopsy.

The left panel of Figure 5 shows the relative magnitudes of the selected fixed
and random effects, while the right panel shows the results of selection for
fixed effects as well as random effects by immune cell sub-class. As described
in Section 4, fixed effects and random effects were considered selection if
their posterior probability of selection was greater than 0.5. Ultimately,
there were 35 covariates for which the fixed effect, random effect, or both
were determined to be significant- see Table S3 of the supplementary materials
for the full list along with immune cell group and effect selection. These
genes came from 6 of the 7 total immune cell subtypes: macrophages, monocytes,
neutrophils, natural killer cells, plasma cells, and T cells. The key results
from the groups with selected genes are detailed in the following paragraphs.

Figure 4: Distribution of spatial interactions within and across biopsies.
Standard deviation of cellular spatial interaction parameter ($\theta$)
plotted against mean $\theta$ parameter at the biopsy level. Also shown are
selected biopsy level plots of $\theta$ by sub-region (darker shades imply
higher interaction) .

#### Neutrophils

The immune cell group with the largest number of effects selected was the
Neutrophil group. Neutrophils play a complex role in the tumor
microenvironment, having been observed to have both pro- and anti-tumorigenic
impacts depending on context (Powell and Huttenlocher, 2016). Notably, one of
the genes for which both the fixed effect and random effect was selected was
PLXNC1, from the Neutrophil gene set. Increased expression of PLXNC1 was
associated with increased tumor-immune interaction between biopsies. This is
consistent with previous biological research that has suggested that PLXNC1
may act as a tumor suppressor, stymying progression and metastasis (Scott et
al., 2009). The combined gene covariate Neutrophil_1 consisted of the genes
S100A8 and S100A9, and was negatively associated with interaction. This is
highly consistent with previous results that indicate that expression of one
or both of these genes is inversely associated with overall survival in
melanoma patients (Pour et al., 2021). We thus find evidence in our results
for the temperamental nature of neutrophil presence in the tumor
microenvironment.

#### T cells

While often treated as a monolith, the T cell family consists of many sub-
types which play vastly differing roles in the TME. For example, while CD8+ T
cells are widely regarded as a key part of the anti-tumorigenic immune
response, regulatory T cells can actually have a pro-tumorigenic impact by
suppressing immune response (Xie et al., 2021).The combined gene covariate T_1
consisted of 72 genes, the full list of which is available in Table S2 of the
supplementary materials, and was positively associated with interaction. While
interpretation is intrinsically more difficult due to the number of genes
being combined, this set most notably contained the genes CD8A (often referred
to as CD8) and CD3E (often referred to as CD3). The presence of one or both of
these genes has been used as a proxy for the presence of cytotoxic lymphocytes
in the context of multiplex imaging (Qin et al., 2022). It is thus
unsurprising that their expression at the biopsy level would be associated
with the presence of such cytotoxic cells, and thus increased interaction
between immune cells and tumor cells. Finally, among the genes where the
random effects were selected was IL23A. IL23A is a pro-inflammatory cytokine
that has been shown to be associated with invasiveness in melanoma (Klein et
al., 2015). While a number of different phenotypes could correspond to
increased heterogeneity within a biopsy, it is intuitive that such
heterogeneity would be associated with increased inflammation and invasiveness
into the surrounding tissue.

Figure 5: Left panel: heatmap of selection probabilities and effect magnitude
for fixed and random effects for selected genes. Posterior mean and selection
probabilities for main effects and random effects. Gene names are colored
according to immune class membership. Right panel: Selection results from gene
expression analysis. Posterior probability of random effect selection plotted
against posterior probability of fixed effect selection. Color indicates class
of immune cells the gene belongs to; shape indicates whether or not the effect
was selected using the median-probability selection model described in Section
4.

#### Monocytes

The role of monocytes in the tumor microenvironment is complex, primarily due
to the fact that monocytes can differentiate into both macrophages and
dendritic cells, each of which has different implications for patient
prognosis. While tumor-associated macrophages are generally associated with
poor patient prognosis, tumor-associated dendritic cells are generally
associated with positive prognosis (Ugel et al., 2021). Amongst the selected
genes in the monocyte group, two stood out based on their prior established
importance in melanoma: HMOX1 and VCAN. HMOX1 has been shown to be genetically
linked to risk of melanoma (Okamoto et al., 2006). Further, downregulating
HMOX1 has been shown to be associated with increased interaction between
natural killer cells and the tumor (Furfaro et al., 2020). This is somewhat
incongruous with our result, since the estimated fixed effect term for HMOX1
is positive; however, since we are essentially adjusting for natural killer
cell presence by proxy (i.e. expression of genes related to natural killer
cells), such an association is not a priori biologically impossible. It has
also been established that VCAN is overexpressed in melanoma relative to
comparable healthy cells and that VCAN is associated with increased
proliferation and in certain melanoma lineages (Touab, 2002; Hernández, 2011).

##  6 Discussion and limitations

In this paper, we propose a method called DreameSpase for integrative analyses
of digital pathology imaging and genomics data. DreameSpase provides a novel
spatially structured regression framework for joint modeling of high
resolution spatial data from non-conformable spaces (i.e. biopsies) and
associated (biopsy-level) covariates (i.e. gene expression data). This is
accomplished by positing a structured regression model with covariate-specific
spatial random effects which are modelled using CAR-based spatial processes.
Furthermore, we induce sparsity by generalizing the traditional spike and slab
priors used for fixed effect selection to the selection of spatial random that
captures the intra-space heterogeneity. We demonstrate via simulation studies
that the model reliably detects inter-space heterogeneity via the fixed effect
selection mechanism, and intra-space heterogeneity via the random effect
selection mechanism. Finally, we demonstrate the utility of DreameSpase in the
analysis of whole slide pathology imaging data by applying it to a set of
images of melanoma biopsies and associated transcriptomics data.

Advances in scale and extent of imaging technologies has led to proliferation
of quantitative imaging data across many diseases, especially cancer. This has
allowed a systematic collation of high-resolution spatial datasets – at scale
– across many patients. These include newer technologies such as spatial
multiplexed tissue imaging that generates high-resolution images to study the
associations of cellular spatial relationships with tumor growth, metastasis,
drug resistance, and patient survival (Van Dam et al., 2022). Our model
engenders joint modeling of such non-conformable spatial data, as well as
assessment of inter and intra-patient heterogeneity and its association with
patient level covariates.

There remain areas of improvement for the model. In our formulation, the
selection of the fixed effects is modeled independently of the selection of
the random effects. This may be a reasonable assumption, but there may be
circumstances in which it may be desirable to allow for the modeling of a
broader class of relationships between the selection of the two types of
effects. Another shortcoming is the formulation of the random effect
selection. The random effect selection suffers from a “piranha problem”
whereby if too many null random effects are included in the model, it can be
challenging for the model to properly select random effects. This is because
each random effect, even if unselected, must have positive variance and thus
account for some small part of the variance of the outcome. In practice, the
estimated values of these random effects tend to be quite small individually,
but large enough in aggregate that if too many are included, the ability of
the model to select true random effects may suffer. In addition to these areas
of improvement, there are also natural potential areas of extension. While we
have limited our inquiry to genomic data, the model could naturally be
extended in the future to multi-omic data (e.g. proteomic, metabolomic data).
While other types of data could simply be included as covariates, it would
perhaps be desirable to introduce group or more structured priors in such a
setting (Zhang et al., 2014). Such priors could be used to select broad types
of covariates (e.g. proteomics vs. genomics) but also potential pathways and
families of genes. We leave these tasks for future exploration.

#### Software and data availability

We have developed a general purpose software package that is available as a
public Github repository. The package includes an efficient implementation of
the Gibbs sampler algorithm discussed in this paper as well as methods to
simulate data from the model. The sampling algorithm is implemented in C++ for
improved efficiency. Data used in the analysis are provided in a separate
Github repository with scripts to reproduce the analysis, tables and figures.
This repository is modular in nature so that individual tables and figures can
be reproduced as necessary.

## References

  * Abousamra et al., (2022) Abousamra, S., Gupta, R., Hou, L., Batiste, R., Zhao, T., Shankar, A., Rao, A., Chen, C., Samaras, D., Kurc, T., and Saltz, J. (2022).  Deep learning-based mapping of tumor infiltrating lymphocytes in whole slide images of 23 types of cancer.  Frontiers in Oncology, 11:806603.
  * Barbieri and Berger, (2004) Barbieri, M. M. and Berger, J. O. (2004).  Optimal predictive model selection.  The Annals of Statistics, 32(3).
  * Baxi et al., (2022) Baxi, V., Edwards, R., Montalto, M., and Saha, S. (2022).  Digital pathology and artificial intelligence in translational medicine and clinical practice.  Modern Pathology, 35(1):23–32.
  * Bowman, (2014) Bowman, F. D. (2014).  Brain imaging analysis.  Annual Review of Statistics and Its Application, 1(1):61–85.
  * Chervoneva et al., (2021) Chervoneva, I., Peck, A. R., Yi, M., Freydin, B., and Rui, H. (2021).  Quantification of spatial tumor heterogeneity in immunohistochemistry staining images.  Bioinformatics, 37(10):1452–1460.
  * Chesney et al., (2018) Chesney, J., Puzanov, I., Collichio, F., Singh, P., Milhem, M. M., Glaspy, J., Hamid, O., Ross, M., Friedlander, P., Garbe, C., Logan, T. F., Hauschild, A., Lebbé, C., Chen, L., Kim, J. J., Gansert, J., Andtbacka, R. H., and Kaufman, H. L. (2018).  Randomized, open-label phase ii study evaluating the efficacy and safety of talimogene laherparepvec in combination with ipilimumab versus ipilimumab alone in patients with advanced, unresectable melanoma.  Journal of Clinical Oncology, 36(17):1658–1667.
  * Clark et al., (2013) Clark, K., Vendt, B., Smith, K., Freymann, J., Kirby, J., Koppel, P., Moore, S., Phillips, S., Maffitt, D., Pringle, M., Tarbox, L., and Prior, F. (2013).  The cancer imaging archive (tcia): Maintaining and operating a public information repository.  Journal of Digital Imaging, 26(6):1045–1057.
  * Falcone et al., (2020) Falcone, I., Conciatori, F., Bazzichetto, C., Ferretti, G., Cognetti, F., Ciuffreda, L., and Milella, M. (2020).  Tumor microenvironment: Implications in melanoma resistance to targeted therapy and immunotherapy.  Cancers, 12(10):2870.
  * Furfaro et al., (2020) Furfaro, A. L., Ottonello, S., Loi, G., Cossu, I., Piras, S., Spagnolo, F., Queirolo, P., Marinari, U. M., Moretta, L., Pronzato, M. A., Mingari, M. C., Pietra, G., and Nitti, M. (2020).  Ho-1 downregulation favors brafv600 melanoma cell death induced by vemurafenib/plx4032 and increases nk recognition.  Int. J. Cancer.
  * George and McCulloch, (1993) George, E. I. and McCulloch, R. E. (1993).  Variable selection via gibbs sampling.  Journal of the American Statistical Association, 88(423):881–889.
  * Goldman et al., (2020) Goldman, M. J., Craft, B., Hastie, M., Repečka, K., McDade, F., Kamath, A., Banerjee, A., Luo, Y., Rogers, D., Brooks, A. N., et al. (2020).  Visualizing and interpreting cancer genomics data via the xena platform.  Nature biotechnology, 38(6):675–678.
  * Hanahan, (2022) Hanahan, D. (2022).  Hallmarks of cancer: New dimensions.  Cancer Discovery, 12(1):31–46.
  * Heindl et al., (2015) Heindl, A., Nawaz, S., and Yuan, Y. (2015).  Mapping spatial heterogeneity in the tumor microenvironment: a new era for digital pathology.  Laboratory Investigation, 95(4):377–384.
  * Hernández, (2011) Hernández, D. (2011).  V3 versican isoform alters the behavior of human melanoma cells by interfering with cd44/erbb-dependent signaling.  Journal of Biological Chemistry, 286(2).
  * Huijben et al., (2023) Huijben, I. A. M., Kool, W., Paulus, M. B., and Van Sloun, R. J. G. (2023).  A review of the gumbel-max trick and its extensions for discrete stochasticity in machine learning.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 45(2).
  * Högmander and Särkkä, (1999) Högmander, H. and Särkkä, A. (1999).  Multitype spatial point patterns with hierarchical interactions.  Biometrics, 55(4):1051–1058.
  * Ibrahim et al., (2011) Ibrahim, J. G., Zhu, H., Garcia, R. I., and Guo, R. (2011).  Fixed and random effects selection in mixed effects models.  Biometrics, 67(2):495–503.
  * Joyner et al., (2020) Joyner, C. N., McMahan, C. S., Tebbs, J. M., and Bilder, C. R. (2020).  From mixed effects modeling to spike and slab variable selection: A bayesian regression model for group testing data.  Biometrics, 76(3):913–923.
  * Kang et al., (2020) Kang, K., Xie, F., Mao, J., Bai, Y., and Wang, X. (2020).  Significance of tumor mutation burden in immune infiltration and prognosis in cutaneous melanoma.  Frontiers in Oncology, 10:573141.
  * Klein et al., (2015) Klein, A., Schwartz, H., Sagi-Assif, O., Meshel, T., Izraely, S., Ben Menachem, S., Bengaiev, R., Ben-Shmuel, A., Nahmias, C., Couraud, P.-O., et al. (2015).  Astrocytes facilitate melanoma brain metastasis via secretion of il-23.  The Journal of pathology, 236(1):116–127.
  * Krishnan et al., (2022) Krishnan, S. N., Mohammed, S., Frankel, T. L., and Rao, A. (2022).  Gawrdenmap: a quantitative framework to study the local variation in cell–cell interactions in pancreatic disease subtypes.  Scientific Reports, 12(1):3708.
  * (22) Li, Q., Wang, X., Liang, F., and Xiao, G. (2019a).  A bayesian mark interaction model for analysis of tumor pathology images.  The Annals of Applied Statistics, 13(3).
  * (23) Li, Q., Wang, X., Liang, F., Yi, F., Xie, Y., Gazdar, A., and Xiao, G. (2019b).  A bayesian hidden potts mixture model for analyzing lung cancer pathology images.  Biostatistics, 20(4):565–581.
  * Li et al., (2021) Li, X., Yang, Y., Huang, Q., Deng, Y., Guo, F., Wang, G., and Liu, M. (2021).  Crosstalk between the tumor microenvironment and cancer cells: A promising predictive biomarker for immune checkpoint inhibitors.  Frontiers in Cell and Developmental Biology, 9:738373.
  * MacNab, (2022) MacNab, Y. C. (2022).  Bayesian disease mapping: Past, present, and future.  Spatial Statistics, 50:100593.
  * Marzagalli et al., (2019) Marzagalli, M., Ebelt, N. D., and Manuel, E. R. (2019).  Unraveling the crosstalk between melanoma and immune cells in the tumor microenvironment.  Seminars in Cancer Biology, 59:236–250.
  * Ni et al., (2019) Ni, Y., Stingo, F. C., and Baladandayuthapani, V. (2019).  Bayesian graphical regression.  Journal of the American Statistical Association, 114(525):184–197.
  * Nirmal et al., (2018) Nirmal, A. J., Regan, T., Shih, B. B., Hume, D. A., Sims, A. H., and Freeman, T. C. (2018).  Immune cell gene signatures for profiling the microenvironment of solid tumors.  Cancer Immunology Research, 6(11):1388–1400.
  * Okamoto et al., (2006) Okamoto, I., Kr, J., Mannhalter, C., Wagner, O., and Pehamberger, H. (2006).  A microsatellite polymorphism in the heme oxygenase-1 gene promoter is associated with risk for melanoma.  International Journal of Cancer.
  * Orozco-Acosta et al., (2022) Orozco-Acosta, E., Adin, A., and Ugarte, M. D. (2022).  Big problems in spatio-temporal disease mapping: methods and software.  Computer Methods and Programs in Biomedicine, 231(arXiv:2201.08323).  arXiv:2201.08323 [stat].
  * Osher et al., (2023) Osher, N., Kang, J., Krishnan, S., Rao, A., and Baladandayuthapani, V. (2023).  Spartin: a bayesian method for the quantification and characterization of cell type interactions in spatial pathology data.  Frontiers in Genetics.
  * Pour et al., (2021) Pour, S. R., de Coaña, Y. P., Demorentin, X. M., Melief, J., Thimma, M., Wolodarski, M., Gomez-Cabrero, D., Hansson, J., Kiessling, R., and Tegner, J. (2021).  Predicting anti-pd-1 responders in malignant melanoma from the frequency of s100a9+ monocytes in the blood.  Journal for immunotherapy of cancer, 9(5).
  * Powell and Huttenlocher, (2016) Powell, D. R. and Huttenlocher, A. (2016).  Neutrophils in the tumor microenvironment.  Trends in Immunology, 37(1):41–52.
  * Qin et al., (2022) Qin, A., Lima, F., Bell, S., Kalemkerian, G. P., Schneider, B. J., Ramnath, N., Lew, M., Krishnan, S., Mohammed, S., Rao, A., and Frankel, T. L. (2022).  Cellular engagement and interaction in the tumor microenvironment predict non-response to pd-1/pd-l1 inhibitors in metastatic non-small cell lung cancer.  Scientific Reports, 12(1):9054.
  * Sadeghi Rad et al., (2021) Sadeghi Rad, H., Monkman, J., Warkiani, M. E., Ladwa, R., O’Byrne, K., Rezaei, N., and Kulasinghe, A. (2021).  Understanding the tumor microenvironment for effective immunotherapy.  Medicinal Research Reviews, 41(3):1474–1498.
  * Saltz et al., (2018) Saltz, J., Gupta, R., Hou, L., Kurc, T., Singh, P., Nguyen, V., Samaras, D., Shroyer, K. R., Zhao, T., Batiste, R., et al. (2018).  Spatial organization and molecular correlation of tumor-infiltrating lymphocytes using deep learning on pathology images.  Cell reports, 23(1):181–193.
  * Scheipl et al., (2012) Scheipl, F., Fahrmeir, L., and Kneib, T. (2012).  Spike-and-slab priors for function selection in structured additive regression models.  Journal of the American Statistical Association, 107(500):1518–1532.
  * Scott et al., (2009) Scott, G. A., McClelland, L. A., Fricke, A. F., and Fender, A. (2009).  Plexin c1, a receptor for semaphorin 7a, inactivates cofilin and is a potential tumor suppressor for melanoma progression.  Journal of Investigative Dermatology, 129(4):954–963.
  * Simiczyjew et al., (2020) Simiczyjew, A., Dratkiewicz, E., Mazurkiewicz, J., Ziȩtek, M., Matkowski, R., and Nowak, D. (2020).  The influence of tumor microenvironment on immune escape of melanoma.  International Journal of Molecular Sciences, 21(21):8359.
  * Sun, (2016) Sun, Y. (2016).  Tumor microenvironment and cancer therapy resistance.  Cancer Letters, 380(1):205–215.
  * Touab, (2002) Touab, M. (2002).  Versican is differentially expressed in human melanoma and may play a role in tumor development.  The American Journal of Pathology, 160(2).
  * Ugel et al., (2021) Ugel, S., Canè, S., De Sanctis, F., and Bronte, V. (2021).  Monocytes in the tumor microenvironment.  Annual Review of Pathology: Mechanisms of Disease, 16(1):93–122.
  * Van Dam et al., (2022) Van Dam, S., Baars, M. J. D., and Vercoulen, Y. (2022).  Multiplex Tissue Imaging: Spatial Revelations in the Tumor Microenvironment.  Cancers, 14(13):3170.
  * Vu et al., (2022) Vu, T., Wrobel, J., Bitler, B. G., Schenk, E. L., Jordan, K. R., and Ghosh, D. (2022).  Spf: A spatial and functional data analytic approach to cell imaging data.  PLOS Computational Biology, 18(6):e1009486.
  * Weinstein et al., (2013) Weinstein, J. N., Collisson, E. A., Mills, G. B., Shaw, K. R. M., Ozenberger, B. A., Ellrott, K., Shmulevich, I., Sander, C., and Stuart, J. M. (2013).  The cancer genome atlas pan-cancer analysis project.  Nature Genetics, 45(10):1113–1120.
  * Xiao et al., (2022) Xiao, X., Guo, Q., Cui, C., Lin, Y., Zhang, L., Ding, X., Li, Q., Wang, M., Yang, W., Kong, Y., and Yu, R. (2022).  Multiplexed imaging mass cytometry reveals distinct tumor-immune microenvironments linked to immunotherapy responses in melanoma.  Communications Medicine, 2(1):131.
  * Xie et al., (2021) Xie, Q., Ding, J., and Chen, Y. (2021).  Role of cd8+ t lymphocyte cells: Interplay with stromal cells in tumor microenvironment.  Acta Pharmaceutica Sinica B, 11(6):1365–1378.
  * Yang et al., (2017) Yang, H., Baladandayuthapani, V., and Morris, J. S. (2017).  Quantile functional regression using quantlets.  arXiv:1711.00031 [stat].  arXiv: 1711.00031.
  * Zhang et al., (2014) Zhang, L., Baladandayuthapani, V., Mallick, B. K., Manyam, G. C., Thompson, P. A., Bondy, M. L., and Do, K.-A. (2014).  Bayesian hierarchical structured variable selection methods with application to molecular inversion probe studies in breast cancer.  Journal of the Royal Statistical Society Series C: Applied Statistics, 63(4):595–620.
  * Zoghi et al., (2023) Zoghi, S., Masoumi, F., and Rezaei, N. (2023).  The immune system.  In Clinical Immunology, pages 1–46. Elsevier.


* [181] K. Xu, Y. Li, R. H. Deng, K. Chen, Deeprefiner: Multi-layer android malware detection system applying deep neural networks, in: 2018 IEEE European Symposium on Security and Privacy (EuroS&P), IEEE, 2018, pp. 473–487.
  * [182] P. Zegzhda, D. Zegzhda, E. Pavlenko, A. Dremov, Detecting android application malicious behaviors based on the analysis of control flows and data flows, in: Proceedings of the 10th International Conference on Security of Information and Networks, 2017, pp. 280–283.
  * [183] Y. Yan, Z. Li, Q. A. Chen, C. Wilson, T. Xu, E. Zhai, Y. Li, Y. Liu, Understanding and detecting overlay-based android malware at market scales, in: Proceedings of the 17th Annual International Conference on Mobile Systems, Applications, and Services, 2019, pp. 168–179.
  * [184] S. Yue, W. Feng, J. Ma, Y. Jiang, X. Tao, C. Xu, J. Lu, Repdroid: an automated tool for android application repackaging detection, in: 2017 IEEE/ACM 25th International Conference on Program Comprehension (ICPC), IEEE, 2017, pp. 132–142.
  * [185] E. Amer, I. Zelinka, S. El-Sappagh, A multi-perspective malware detection approach through behavioral fusion of api call sequence, Computers & Security 110 (2021) 102449.
  * [186] J. Feng, L. Shen, Z. Chen, Y. Wang, H. Li, A two-layer deep learning method for android malware detection using network traffic, IEEE Access 8 (2020) 125786–125796.
  * [187] H. Wang, W. Zhang, H. He, You are what the permissions told me! android malware detection based on hybrid tactics, Journal of Information Security and Applications 66 (2022) 103159.
  * [188] S. Bhandari, R. Panihar, S. Naval, V. Laxmi, A. Zemmari, M. S. Gaur, Sword: semantic aware android malware detector, Journal of information security and applications 42 (2018) 46–56.
  * [189] L. Massarelli, L. Aniello, C. Ciccotelli, L. Querzoni, D. Ucci, R. Baldoni, Android malware family classification based on resource consumption over time, in: 2017 12th International Conference on Malicious and Unwanted Software (MALWARE), IEEE, 2017, pp. 31–38.
  * [190] Y. Zhao, G. Xu, Y. Zhang, Hfa-md: An efficient hybrid features analysis based android malware detection method, in: International Conference on Heterogeneous Networking for Quality, Reliability, Security and Robustness, Springer, 2017, pp. 248–257.
  * [191] M. Alhanahnah, Q. Yan, H. Bagheri, H. Zhou, Y. Tsutano, W. Srisa-An, X. Luo, Dina: Detecting hidden android inter-app communication in dynamic loaded code, IEEE Transactions on Information Forensics and Security 15 (2020) 2782–2797.
  * [192] P. Xu, Y. Zhang, C. Eckert, A. Zarras, Hawkeye: cross-platform malware detection with representation learning on graphs, in: International Conference on Artificial Neural Networks, Springer, 2021, pp. 127–138.
  * [193] F. Shahzad, S. Bhatti, M. Shahzad, M. Farooq, In-execution malware detection using task structures of linux processes, in: 2011 IEEE International Conference on Communications (ICC), 2011, pp. 1–6.
  * [194] A. Shalaginov, L. Øverlier, A novel study on multinomial classification of x86/x64 linux elf malware types and families through deep neural networks, in: Malware Analysis Using Artificial Intelligence and Deep Learning, Springer, 2021, pp. 437–453.
  * [195] E. Cozzi, M. Graziano, Y. Fratantonio, D. Balzarotti, Understanding linux malware, in: 2018 IEEE Symposium on Security and Privacy (SP), 2018, pp. 161–175.
  * [196] J. Carrillo-Mondéjar, J. L. Martínez, G. Suarez-Tangil, Characterizing linux-based malware: Findings and recent trends, Future Generation Computer Systems 110 (2020) 267–281.
  * [197] K. Monnappa, Automating linux malware analysis using limon sandbox, Black Hat Europe 2015 (2015) IV–A.
  * [198] A. Yazdinejad, H. HaddadPajouh, A. Dehghantanha, R. M. Parizi, G. Srivastava, M.-Y. Chen, Cryptocurrency malware hunting: A deep recurrent neural network approach, Applied Soft Computing 96 (2020) 106630.
  * [199] E. Snow, M. Alam, A. Glandon, K. Iftekharuddin, End-to-end multimodel deep learning for malware classification, in: 2020 International Joint Conference on Neural Networks (IJCNN), 2020, pp. 1–7.
  * [200] M. Khan, D. Baig, U. S. Khan, A. Karim, Malware classification framework using convolutional neural network, in: 2020 International Conference on Cyber Warfare and Security (ICCWS), 2020, pp. 1–7.
  * [201] S. K. Smmarwar, G. P. Gupta, S. Kumar, Design of a fused triple convolutional neural network for malware detection: A visual classification approach, in: International Conference on Advances in Computing and Data Sciences, Springer, 2021, pp. 279–289.
  * [202] F. O. Catak, J. Ahmed, K. Sahinbas, Z. H. Khand, Data augmentation based malware detection using convolutional neural networks, PeerJ Computer Science 7 (2021) e346.
  * [203] A. Bensaoud, J. Kalita, Deep multi-task learning for malware image classification, Journal of Information Security and Applications 64 (2022) 103057\.
  * [204] Y. Ding, M. Shao, C. Nie, K. Fu, An efficient method for generating adversarial malware samples, Electronics 11 (1) (2022) 154.
  * [205] K. Marsh, H. Haddadpajouh, Ransomware threat detection: A deep learning approach, in: Handbook of Big Data Analytics and Forensics, Springer, 2022, pp. 253–269.
  * [206] A. Tekerek, M. M. Yapici, A novel malware classification and augmentation model based on convolutional neural network, Computers & Security 112 (2022) 102515\.
  * [207] I. Obaidat, M. Sridhar, K. M. Pham, P. H. Phung, Jadeite: A novel image-behavior-based approach for java malware detection using deep learning, Computers & Security 113 (2022) 102547.
  * [208] I. Ben Abdel Ouahab, L. Elaachak, M. Bouhorma, Image-based malware classification using multi-layer perceptron, in: Networking, Intelligent Systems and Security, Springer, 2022, pp. 453–464.
  * [209] F. O. Catak, A. F. Yazı, O. Elezaj, J. Ahmed, Deep learning based sequential model for malware analysis using windows exe api calls, PeerJ Computer Science 6 (2020) e285.
  * [210] M. Ring, D. Schlör, S. Wunderlich, D. Landes, A. Hotho, Malware detection on windows audit logs using lstms, Computers & Security 109 (2021) 102389.
  * [211] W. R. Aditya, Girinoto, R. B. Hadiprakoso, A. Waluyo, Deep learning for malware classification platform using windows api call sequence, in: 2021 International Conference on Informatics, Multimedia, Cyber and Information System (ICIMCIS, 2021, pp. 25–29.
  * [212] S. K. J. Rizvi, W. Aslam, M. Shahzad, S. Saleem, M. M. Fraz, Proud-mal: static analysis-based progressive framework for deep unsupervised malware classification of windows portable executable, Complex & Intelligent Systems 8 (1) (2022) 673–685.
  * [213] J.-Y. Kim, S.-B. Cho, Obfuscated malware detection using deep generative model based on global/local features, Computers & Security 112 (2022) 102501.
  * [214] L. Xin, L. Chao, L. He, Malicious code detection method based on image segmentation and deep residual network resnet, in: 2021 International Conference on Computer Engineering and Application (ICCEA), 2021, pp. 473–480.
  * [215] R. Damaševičius, A. Venčkauskas, J. Toldinas, Š. Grigaliūnas, Ensemble-based classification using neural networks and machine learning models for windows pe malware detection, Electronics 10 (4) (2021) 485.
  * [216] Y. Jian, H. Kuang, C. Ren, Z. Ma, H. Wang, A novel framework for image-based malware detection with a deep neural network, Computers & Security 109 (2021) 102400.
  * [217] A. Mallik, A. Khetarpal, S. Kumar, Conrec: malware classification using convolutional recurrence, Journal of Computer Virology and Hacking Techniques (2022) 1–17.
  * [218] D. Mane, P. Kumbharkar, S. B. Javheri, R. Moorthy, An adaptable ensemble architecture for malware detection, in: International Conference on Innovative Computing and Communications, Springer, 2022, pp. 647–659.
  * [219] L. Li, Y. Ding, B. Li, M. Qiao, B. Ye, Malware classification based on double byte feature encoding, Alexandria Engineering Journal 61 (1) (2022) 91–99.
  * [220] Y. Fang, Y. Zeng, B. Li, L. Liu, L. Zhang, Deepdetectnet vs rlattacknet: An adversarial method to improve deep learning-based static malware detection model, Plos one 15 (4) (2020) e0231626.
  * [221] I. Alodat, M. Alodat, Detection of image malware steganography using deep transfer learning model, in: Proceedings of International Conference on Data Science and Applications, Springer, 2022, pp. 323–333.
  * [222] D. Carlin, A. Cowan, P. O’Kane, S. Sezer, The effects of traditional anti-virus labels on malware detection using dynamic runtime opcodes, IEEE Access 5 (2017) 17742–17752.
  * [223] A. Pektaş, T. Acarman, Learning to detect android malware via opcode sequences, Neurocomputing 396 (2020) 599–608.
  * [224] Z. Ma, H. Ge, Z. Wang, Y. Liu, X. Liu, Droidetec: Android malware detection and malicious code localization through deep learning, arXiv preprint arXiv:2002.03594.
  * [225] A. Lakshmanarao, M. Shashi, Android malware detection with deep learning using rnn from opcode sequences., International Journal of Interactive Mobile Technologies 16 (1).
  * [226] R. Feng, J. Q. Lim, S. Chen, S.-W. Lin, Y. Liu, Seqmobile: An efficient sequence-based malware detection system using rnn on mobile devices, in: 2020 25th International Conference on Engineering of Complex Computer Systems (ICECCS), IEEE, 2020, pp. 63–72.
  * [227] A. Darwaish, F. Naït-Abdesselam, Rgb-based android malware detection and classification using convolutional neural network, in: GLOBECOM 2020-2020 IEEE Global Communications Conference, IEEE, 2020, pp. 1–6.
  * [228] R. Feng, S. Chen, X. Xie, G. Meng, S.-W. Lin, Y. Liu, A performance-sensitive malware detection system using deep learning on mobile devices, IEEE Transactions on Information Forensics and Security 16 (2020) 1563–1578.
  * [229] G. Iadarola, F. Martinelli, F. Mercaldo, A. Santone, Towards an interpretable deep learning model for mobile malware detection and family identification, Computers & Security 105 (2021) 102198.
  * [230] J. Kim, Y. Ban, E. Ko, H. Cho, J. H. Yi, Mapas: a practical deep learning-based android malware detection system, International Journal of Information Security (2022) 1–14.
  * [231] K. Kong, Z. Zhang, Z.-Y. Yang, Z. Zhang, Fcscnn: Feature centralized siamese cnn-based android malware identification, Computers & Security 112 (2022) 102514\.
  * [232] I. Almomani, A. Alkhayer, W. El-Shafai, An automated vision-based deep learning model for efficient detection of android malware attacks, IEEE Access 10 (2022) 2700–2720.
  * [233] H.-I. Kim, M. Kang, S.-J. Cho, S.-I. Choi, Efficient deep learning network with multi-streams for android malware family classification, IEEE Access 10 (2022) 5518–5532.
  * [234] Y. Ban, S. Lee, D. Song, H. Cho, J. H. Yi, Fam: Featuring android malware for deep learning-based familial analysis, IEEE Access 10 (2022) 20008–20018.
  * [235] C. Yuan, J. Cai, D. Tian, R. Ma, X. Jia, W. Liu, Towards time evolved malware identification using two-head neural network, Journal of Information Security and Applications 65 (2022) 103098.
  * [236] S. Garg, N. Baliyan, M2vmapper: Malware-to-vulnerability mapping for android using text processing, Expert Systems with Applications 191 (2022) 116360.
  * [237] M. Almahmoud, D. Alzu’bi, Q. Yaseen, Redroiddet: Android malware detection based on recurrent neural network, Procedia Computer Science 184 (2021) 841–846.
  * [238] S. K. Sasidharan, C. Thomas, Memdroid-lstm based malware detection framework for android devices, in: 2021 IEEE Pune Section International Conference (PuneCon), 2021, pp. 1–6.
  * [239] S. Sabhadiya, J. Barad, J. Gheewala, Android malware detection using deep learning, in: 2019 3rd International Conference on Trends in Electronics and Informatics (ICOEI), IEEE, 2019, pp. 1254–1260.
  * [240] P. Feng, J. Ma, T. Li, X. Ma, N. Xi, D. Lu, Android malware detection via graph representation learning, Mobile Information Systems 2021.
  * [241] Y. Wu, J. Shi, P. Wang, D. Zeng, C. Sun, Deepcatra: Learning flow-and graph-based behaviors for android malware detection, arXiv preprint arXiv:2201.12876.
  * [242] S. Mahdavifar, D. Alhadidi, A. Ghorbani, et al., Effective and efficient hybrid android malware classification using pseudo-label stacked auto-encoder, Journal of Network and Systems Management 30 (1) (2022) 1–34.
  * [243] H. Gao, S. Cheng, W. Zhang, Gdroid: Android malware detection and classification with graph convolutional network, Computers & Security 106 (2021) 102264.
  * [244] T. Lu, Y. Du, L. Ouyang, Q. Chen, X. Wang, Android malware detection based on a hybrid deep learning model, Security and Communication Networks 2020.
  * [245] E. Mariconti, L. Onwuzurike, P. Andriotis, E. De Cristofaro, G. Ross, G. Stringhini, Mamadroid: Detecting android malware by building markov chains of behavioral models, arXiv preprint arXiv:1612.04433.
  * [246] S. Laudanna, C. A. Visaggio, et al., Gang-mam: Gan based engine for modifying android malware, arXiv preprint arXiv:2109.13297.
  * [247] S. Seraj, S. Khodambashi, M. Pavlidis, N. Polatidis, Hamdroid: permission-based harmful android anti-malware detection using neural networks, Neural Computing and Applications (2022) 1–10.
  * [248] D. Palmer, Linux malware attacks are on the rise, and businesses aren’t ready for it — zdnet, https://www.zdnet.com/article/linux-malware-attacks-are-on-the-rise-and-businesses-arent-ready-for-it/, (Accessed on 05/22/2022) (2022).
  * [249] VMware, Threat report – exposing malware in linux-based multi-cloud environments - vmware security blog - vmware, https://blogs.vmware.com/security/2022/02/2022-vmware-threat-report-exposing-malware-in-linux-based-multi-cloud-environments.html, (Accessed on 05/22/2022) (2022).
  * [250] T. Landman, N. Nissim, Deep-hook: A trusted deep learning-based framework for unknown malware detection and classification in linux cloud environments, Neural Networks 144 (2021) 648–685.
  * [251] T. Panker, N. Nissim, Leveraging malicious behavior traces from volatile memory using machine learning methods for trusted unknown malware detection in linux cloud environments, Knowledge-Based Systems 226 (2021) 107095.
  * [252] C. Hwang, J. Hwang, J. Kwak, T. Lee, Platform-independent malware analysis applicable to windows and linux environments, Electronics 9 (5) (2020) 793.
  * [253] R. C., Linux malware on a rise reaching all-time high in h1 2022 - atlas vpn, https://atlasvpn.com/blog/linux-malware-on-a-rise-reaching-all-time-high-in-h1-2022, (Accessed on 08/18/2022) (2022).
  * [254] C. Zhang, S. Bengio, M. Hardt, B. Recht, O. Vinyals, Understanding deep learning (still) requires rethinking generalization, Communications of the ACM 64 (3) (2021) 107–115.
  * [255] D. Patel, P. Sastry, Memorization in deep neural networks: Does the loss function matter?, in: Pacific-Asia Conference on Knowledge Discovery and Data Mining, Springer, 2021, pp. 131–142.
  * [256] Appbrain, Android and google play statistics, development resources and intelligence — appbrain, https://www.appbrain.com/stats, (Accessed on 05/19/2022) (2022).
  * [257] BuildFire, Mobile app download statistics & usage statistics (2022) - buildfire, https://buildfire.com/app-statistics/, (Accessed on 05/19/2022) (2022).
  * [258] M. D. F. De Grazia, I. Stoianov, M. Zorzi, Parallelization of deep networks., in: ESANN, Citeseer, 2012.
  * [259] R. Mu, X. Zeng, A review of deep learning research, KSII Transactions on Internet and Information Systems (TIIS) 13 (4) (2019) 1738–1764.
  * [260] O. Suciu, S. E. Coull, J. Johns, Exploring adversarial examples in malware detection, in: 2019 IEEE Security and Privacy Workshops (SPW), IEEE, 2019, pp. 8–14.
  * [261] K. Grosse, N. Papernot, P. Manoharan, M. Backes, P. McDaniel, Adversarial examples for malware detection, in: European symposium on research in computer security, Springer, 2017, pp. 62–79.
  * [262] F. Kreuk, A. Barak, S. Aviv-Reuven, M. Baruch, B. Pinkas, J. Keshet, Deceiving end-to-end deep learning malware detectors using adversarial examples, arXiv preprint arXiv:1802.04528.
  * [263] H. S. Anderson, A. Kharkar, B. Filar, D. Evans, P. Roth, Learning to evade static pe machine learning malware models via reinforcement learning, arXiv preprint arXiv:1801.08917.
  * [264] D. Park, H. Khan, B. Yener, Generation & evaluation of adversarial examples for malware obfuscation, in: 2019 18th IEEE International Conference On Machine Learning And Applications (ICMLA), IEEE, 2019, pp. 1283–1290.
  * [265] D. Manky, How linux became the new bullseye for bad guys — securityweek, https://www.securityweek.com/how-linux-became-new-bullseye-bad-guys, (Accessed on 05/22/2022) (2022).
  * [266] R. Lakshmanan, Microsoft warns rise in xorddos malware targeting linux devices, https://thehackernews.com/2022/05/microsoft-warns-rise-in-xorddos-malware.html, (Accessed on 05/24/2022) (2022).