# Near-Field Radiation Exposure Control in Slot-Loaded Microstrip Antenna: A
Characteristic Mode Approach

Sandip Ghosal, Arijit De, Raed M. Shubair, and Ajay Chakrabarty S. Ghosal, A.
De and A. Chakrabarty are with the Department of Electronics and Electrical
Communication Engineering, Indian Institute of Technology, Kharagpur, West
Bengal 721302, India (e-mail: sgmw@iitkgp.ac.in).Raed M. Shubair is with
Research Laboratory of Electronics, Massachusetts Institute of Technology
(MIT), Cambridge, MA 02139 USA (e-mail: rshubair@mit.edu).

###### Abstract

Microstip antenna topology is commonly loaded with a narrow slot to manipulate
the resonance frequency or impedance bandwidth. However, the tuning of the
resonance frequency or impedance bandwidth results in the variation of the
current and field distributions. In this regard, this work adopts the concept
of characteristic modes to gain an initial understanding of the perturbation
mechanism of the rectangular patch when loaded with a slot. The performance of
microstrip antennas with finite ground plane is then studied using full-wave
simulation. It has been found that the distribution of the induced current
density is highly dependent on the orientation of the slot The incorporation
of a narrow slot suppresses the nearby orthogonal eigen mode and, as a
consequence, the radiation behaviour is affected. Specifically, in the
presence of biological tissues in the near-field region, both antenna input
impedance properties and the realized gain are dependent on the slot
orientation. Different examples are included for understanding the impact of
slot loading on the energy absorption by biological tissues, by calculating
the the specific absorption rate (SAR). The proposed analysis facilitates the
design of miniaturized antenna geometries for biomedical applications via
systematic loading of narrow slots.

##  I Introduction

The last two decades have witnessed an exponential growth and tremendous
developments in antenna technologies and techniques, and their associated
applications, such as those reported in [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44] . Being an integral
part of the cellular phone, antenna’s near-field exposure on biological
tissues is a major concern for the antenna designers. RF exposure induced by
the GSM mobile phone has been studied using statistical modelling in [45]. A
numerical technique has also been proposed to determine the RF exposure
effects on a human head due to a cellphone equipped with a dual-band monopole-
helix antenna [46]. In a similar way, various statistical and experimental
studies were carried out to assess the impact of mobile radiation on nearby
biological tissues [47] [48] [49]. Recently, modal analysis technique has been
adopted for predicting the frequency of maximum RF exposure [50] [51].
Therefore, the antenna geometry should be designed in a way such that it
maintains a certain threshold of RF energy absorption by the biological
tissues commonly characterized by the SAR value [52] [53].

Antenna design in the presence of nearby biological tissues involves three
main aspects. The first is to select the optimum frequency for maximum power
transfer, as rigorously studied in [54]. The second is to focus the radiated
field or SAR distribution towards a desired direction [55] [56]. The third and
last aspect, is to miniaturize the design while maintaining the desired
radiation characteristics in the desired frequency of operation.
Miniaturization of the design is achieved by incorporating electromagnetic
band gap (EBG) structures [57] or slot elements [58] in the topology of the
antenna. The latter technique is the focus of this report.

Slot loading of is a popular technique used by the antenna designers to
manipulate resonance frequency [59], enhance impedance bandwidth [60], or
reduce side-lobe level [61]. The principle behind slot loading technique is to
manipulate the radiation behaviour which is directly controlled by the current
density vector induced on the surface or volume of the radiating antenna. In
this regard, a proper analysis is required to study the effect of slot loading
on the radiated field and SAR. In the literature, slot loading phenomena has
been discussed following three approaches: cavity modelling [61], method of
moment (MoM) technique [62], and transmission line theory [58]. Recently, a
modal approach based on characterisitc mode analysis (CMA) has been employed
to study the effect of slot loading [63] [64]. Characteristic mode theory was
initially developed to determine the induced current vector using generalized
eigen-decomposition of the MoM impedance matrix. If Galerkin’s method is
followed in the electric field integral equation (EFIE) formulation, the MoM
impedance matrix becomes symmetric [65]. Generalized characteristic equation
of the MoM matrix provides the eigen-values and the characteristic current
eigen-modes. As shown in [65] and [51], these eigen-values can predict the
resonance frequency of the antenna structure. The amount of outward radiated
power is also directly dependent on the eigen-modes and eigen-values. Unlike
the conventional MoM technique, computation of these eigen parameters is
independent of the excitation vector. It depends only on the operating
frequency, and on the physical and geometrical properties of the radiating
structure. Hence, with the knowledge of excitation-independent eigen-current
distribution, one can perturb the radiating structure topology to achieve a
desired radiation pattern. Such advantageous scopes of the CMA are found to be
recently exploited in controlling the higher order modes [66] and radar cross
section [67] respectively. This report aims to extend the previously reported
modal approach of [63] [64], for perturbing the horizontal radiation in the
near field of the radiating structure, such that the RF exposure on the on the
nearby biological tissues is adjusted.

The report is organized as follows: Section II starts by describing the
characteristic mode analysis, then considers a rectangular geometry to study
the perturbation of modal current density due to the loading of narrow slots.
Section III explores the method further with full-wave examples and
corresponding analysis of the SAR distribution Finally, conclusions are given
in Section IV..

##  II Modal Analysis

For the sake of initial understanding, let us consider a perfect electrically
conducting (PEC) radiating element with no internal resonance. The outward
unit vector normal to the surface of the antenna is denoted by $\hat{n}$.
$\text{E}^{i}$ is the excitation field vector of the antenna. J is the surface
current density vector induced on the surface $S$ of the antenna. Following
the field boundary conditions on the surface $S$ [65],

$\displaystyle\hat{n}\times\text{E}^{i}$ $\displaystyle=-\hat{n}\times[j\omega
A(\text{J}(r^{\prime})+\nabla\phi(\text{J}(r^{\prime})]$ (1a) $\displaystyle
A(\text{J}(r^{\prime})$
$\displaystyle=\mu\oint_{S}\text{J}(r^{\prime})G(r,r^{\prime})\,ds$ (1b)
$\displaystyle\phi(\text{J}(r^{\prime})$
$\displaystyle=\frac{-1}{j\omega\epsilon}\oint_{S}\nabla^{\prime}\text{J}(r^{\prime})G(r,r^{\prime})\,ds$
(1c) $\displaystyle G(r,r^{\prime})$
$\displaystyle=\frac{e^{-jk|r-r^{\prime}|}}{4\pi|r-r^{\prime}|}$ (1d)

The permittivity, permeability and wave-number of the free-space medium are
represented by $\mu$, $\epsilon$ and k respectively. The integral equation of
(1a)-(1d) can be converted into equivalent matrix form using the method of
moment (MoM) technique [65] as,

$[Z]_{M\times M}[I]_{M\times 1}=[V]_{M\times 1}$ (2)

The MoM impedance matrix is denoted by $[Z]$ and the column vector $[V]$
consists of the weighting coefficients generated due to the inner product of
the excitation vector $\text{E}^{i}$ and $M$ number of basis functions
$\text{f}_{m}$. Using the CMA of [65], the weighting vector $[I]$ of the
induced current can be written as

$\displaystyle[I]$
$\displaystyle=\sum_{m=1}^{M}\alpha_{m}[I_{m}],~{}\text{where}~{}\alpha_{m}=\frac{[I_{m}]^{H}[V]}{(1+j\lambda_{m})[I_{m}]^{H}[R][I_{m}]}$
(3)

and the $m^{th}$ characteristic current mode $\text{J}_{m}$ will be

$\displaystyle\text{J}_{m}$
$\displaystyle=\sum_{i=1}^{M}I_{m}(i)\text{f}_{i}(r^{\prime}_{m}),r^{\prime}_{m}\in
S$ (4)

The $m^{th}$ eigenvector $[I_{m}]$ and eigenvalue $\lambda_{m}$ can be
obtained through the generalised eigen decomposition of the impedance matrix
$[Z]$ as,

$\displaystyle[X][I_{m}]=\lambda_{m}[R][I_{m}],~{}\text{where}~{}[Z]=[R]+j[X]$
(5)

If Galerkin’s matching is followed in computing the $[Z]$ matrix, both $[R]$
and $[X]$ will be real symmetric. So the eigen pair $([I_{m}],\lambda_{m})$ of
(5) will be necessarily real valued. However it can be noted from [51] that
the lossy structures will have complex eigenvectors and eigenvalues. Following
the orthogonal properties of the eigenvectors of (5), it can be written that

$\displaystyle[I_{m}]^{H}[R][I_{n}]=\delta_{mn}[I_{m}]^{H}[R][I_{m}]$ (6a)
$\displaystyle[I_{m}]^{H}[X][I_{n}]=\lambda_{m}\delta_{mn}[I_{m}]^{H}[R][I_{m}]$
(6b) $\displaystyle\delta_{mn}=1~{}\text{if}~{}m=n,~{}\text{else}~{}0$ (6c)

The Hermitian transpose operation is denoted by “H”. for the lossless
structures, it will turn out to be as transposition operation. The modal
significance of the $m^{th}$ eigenmode is denoted by

$MS_{m}=\frac{1}{\sqrt{1+\lambda_{m}^{2}}}.$ (7)

Following the complex -power balance relation of [65], the power radiated by
the $m^{th}$ mode can be defined as,

$P_{R,m}=|\alpha_{m}|^{2}=MS_{m}(\frac{[I_{m}]^{H}[V]}{[I_{m}]^{H}[R][I_{m}]})^{2}$
(8)

Total radiated power at a certain frequency,say,$f_{0}$ will be a weighted sum
of all individual modal power,

$P_{R}=\sum_{m=1}^{M}P_{R,m}=\sum_{m=1}^{M}|\alpha_{m}|^{2}$ (9)

Now, let’s come back to the first two design objectives as stated in the
initial part of the previous section. It can be noted from (8) that the
outward radiated power $P_{R,m}$ is inversely related to the eigenvalue
$\lambda_{m}$ of the corresponding mode. As $\lambda_{m}$ tends towards 0,
corresponding modal significance $MS_{m}$ of (7) and radiated power $P_{R,m}$
of (9) will be higher. Similarly when the antenna is used in the receiving
end, real power of the pointing vectors provides $P_{R,m}$ as the received
real power. So, if the excitation vector is restricted to a very finite region
compared to the total geometry of the antenna structure, the optimal frequency
of maximum radiation can be approximated following the minima of the
corresponding eigenvalue.

Figure 1: (a) Reference plate with $L=50$ mm, $W=40$ mm. (b)The plate with
narrow slot. For X-directed slot, $L_{s}=30$ mm and $W_{s}=0.2$ mm. For
Y-directed slot, $L_{s}=0.2$ mm and $W_{s}=30$ mm.

On the other hand, total radiated power consists of the contribution from all
modes in (9). To align the field distribution along a certain direction or to
achieve a desired SAR distribution, it is required to suppress the adjacent
undesired modes. Let the $i^{th}$ mode generates the desired beam pattern at
the resonating frequency of $\omega_{i}$. Henceforth the desired radiated
modal power is represented by $|\alpha_{m}|^{2}$ where $m=i$. Alternately, the
undesired component of the radiated power are contributed by the other
additional modes, $\sum_{m=1}^{M}|\alpha_{m}|^{2}$ for $m\neq i$. So the
control over undesired radiation needs to start from the search of the
characteristic mode(/s) whose modal radiated power is comparable to the
desired mode’s power. Physically the eigenvalue parameter $\lambda_{i}$
denotes the proximity towards resonance of the $i^{th}$ characteristic mode.
Ideal resonance criterion is defined by zero reactive power,in other words,
$\lambda_{m}=0$. So, the dominant eigenmodes are commonly determined using the
modal significance parameter $MS_{m}$ of (7). The next discussion shows how
the modal characteristics vary due to loading of narrow slot. For the purpose
of illustration, let us consider a rectangular PEC plate in Fig.1. The plate
loaded with narrow slot is shown in Fig.1.

(a) Mode 1, Ref.

(b) Mode 2, Ref.

(c) Mode 1, X-slot

(d) Mode 2, X-slot

(e) Mode 1, Y-slot

(f) Mode 2, Y-slot

Figure 2: Normalized distribution of the $1^{st}$ and $2^{nd}$ eigencurrent
modes at $f=2.4$ GHz.

The eigenmodes have been computed in MATLAB [68] following [65] and [69] with
RWG basis functions. Sorting the modes following (7), the $1^{st}$ two
dominant eigen current modes at $f=2.4$ GHz have been plotted in Figs.2(a) and
2(b). Keeping both the current distribution in same scale, the eigencurrent
modes have been normalised with respect to its maximum amplitude. This leads
to the presence of the additional term $[I_{m}]^{H}[R][I_{m}]$ in defining
$\alpha_{m}$ of (3). As shown in Figs.2(a) and 2(b), the $1^{st}$ and $2^{nd}$
eigenmodes are found to be linearly polarized along X and Y-direction
respectively.
In the next stage, narrow slot has been incorporated with its wider axis along
the X-direction where $Ws=0.2$ mm and $Ls=30$ mm. The perturbed current
distribution of the first two dominant modes at $f=2.4$ GHz are shown in Figs.
2(c) and 2(d). Major contribution of the slot is reflected in the second mode
of Fig. 2(d). It can be noted from the comparison between Fig.2(b) and 2(d)
that the X-directed slot converts the nearest Y-polarized edge mode of
Fig.2(b) into a slot mode of Fig.2(d) where the current distribution hold
significant magnitude only along two narrow ends of the slot. The X-directed
slot has minimal impact on the first mode in Fig.2(c) whose orientation is
parallel to the wider axis of the slot. For studying the complementary case,
the slot has been aligned along the Y-direction and the resultant modes have
been shown in Figs. 2(e) and 2(f). Similar to the earlier case, the
Y-polarised edge mode of Fig.2(a) is transformed into the Y-polarised edge
mode in Fig.2(e). Minimal variation is noted in the X-polarised edge mode of
Fig.2(f).

Figure 3: Field observation plane at $z=5$ mm height from the radiating plate.

(a) $20log_{10}|\vec{E}_{1}|$ along X-axis

(b) $20log_{10}|\vec{E}_{1}|$ along Y-axis

(c) $20log_{10}|\vec{E}_{2}|$ along X-axis

(d) $20log_{10}|\vec{E}_{2}|$ along X-axis

Figure 4: Near field distribution due to the first two dominant characteristic
current modes at $f=2.4$ GHz.

As the loading of slot modify the current distribution, corresponding field
distribution will also change. The SAR distribution on a nearby object is
directly related to the electric field as

$\displaystyle\text{SAR}=\frac{\sigma}{\rho}|\text{E}|^{2}$ (10)

The nearby object will experience the incident field E contributed by the
antenna structure. For example while holding a cell phone like Fig. 3, the
horizontal radiation will influence the biological tissues. To understand the
impacts on the near field variation along the horizontal direction, let us
consider a line of observation at an offset distance from the edges of the PEC
plate (i.e., the radiating source). Using infinitesimal dipole model of [69],
the electric field $20log_{10}|\text{E}|$ has been calculated for the
corresponding eigencurrent modes. It can be noted from Figs. 4(a) and 4(b)
that the X-directed slot enhances the X-polarised $1^{st}$ mode’s near field
density where as the Y-directed slot reduces the corresponding modal strength.
Similarly, the Y-polarised mode is found to be weakened by the X-directed
narrow slot in Figs. 4(c) and 4(d). Such variation in the near field behaviour
due to loading of narrow slot will influence the SAR distribution of the
neighbouring biological tissue following (10). Therefore, the radiating
element can be suitably designed based the characteristic mode analysis to
minimize the undesired radiation. Till now, the modal analysis does not
consider the information of the excitation vector. The next section uses the
same radiating elements of Figs. 1 and 1 with finite microstip ground plane to
verify the suitability of the slot loading method using full-wave simulator
Ansys HFSS [70].

Figure 5: Antenna with the adjacent tissue layer.

##  III Full-wave Analysis

For the full-wave study, two different case studies have been considered with
probe-feeding and microstrip line feeding arrangement. For each case, the
patch dimensions have been assumed to be same as of Figs. 1 and 1. The patch
surface is lying 1.57 mm above the lossless Roger substrate with the
dielectric permittivity of 2.2. The tissue is placed parallel to the XZ-plane
with the Y-separation of 2 mm from the substrate edge of the antenna as shown
in Fig. 5.

###  III-A Probe-fed Patch

The dimension of the ground plane is assumed to be 83.6 mm $\times$ 71.5 mm
for tuning the reference antenna around 2.4 GHz using Ansys Design Kit. The
antenna has been excited by co-axial probe feeding at ($x_{f}=2$ mm,
$y_{f}=8.05$ mm, $z_{f}=1.57$ mm). The current distribution induced on the
patch surface of the reference antenna is shown in Fig. 6 where the the
current is found to be Y-polarised. Following the previous modal discussion,
if it is desired to suppress this Y-polarised current density, the narrow slot
needs to be aligned along the orthogonal direction, i.e., along the X-axis. So
for the sake of comparison, two examples have been studied with the X-directed
and Y-directed slots on the patch surface. The feed location has been taken as
same of the reference probe-fed patch. As obvious in Fig. 6, the orthogonal
orientation of the slot can perturb the Y-polarised edge mode of Fig. 6 into
the slot mode. On the other hand, as the slot is oriented parallel to the
Y-polarised edge mode, it has negligible impact on the current distribution in
Fig. 6. As the slot perturbs the induced current density vector, the impedance
characteristics also change in Fig. 6. The change in the equivalent circuit
parameters can be computed using the existing analysis of [58] and [63]. The
next discussion shows how the perturbation of the current vector becomes
reflected in the near field and and corresponding SAR distribution.

Figure 6: Normalized distribution of the characteristic current modes at
$f=2.4$ GHz. (c) $1^{st}$ and (d) $2^{nd}$ modes of the plate with X-directed
slot.

(a) $20log_{10}|\vec{E}_{1}|$ along X-axis

(b) $20log_{10}|\vec{E}_{1}|$ along Y-axis

Figure 7: Near field distribution due to the probe-fed patch at $f=2.4$ GHz.

(a) $\phi=0^{\circ}$-plane

(b) $\phi=0^{\circ}$-plane

Figure 8: Far field distribution due to the probe-fed patch at $f=2.4$ GHz.

To understand the change in horizontal radiation, the near field at the same
offset distance from two edges are reported in Figs. 7(a) and 7(b). Similar to
the previous observation of Figs. 4(a)–4(d), the near field strength along the
horizontal directed is found to be reduced due to loading of the narrow slot
of orthogonal orientation. In the next stage, the tissue layer is placed
parallel to the Y-edge of the antenna as shown in Fig. 5. The separation
distance between the antenna and the tissue layer along the Y-axis is assumed
to be 2 mm.

TABLE I: Specification of the tissue layers

Tissue | $\mathbf{\sigma(S/m)}$ | $\mathbf{\epsilon_{r}}$ | $\mathbf{\rho(kg/m^{3})}$ | d(mm)
---|---|---|---|---
skin | $1.43$ | $38.1$ | 1100 | 1.5
fat | $0.1$ | $5.29$ | 916 | 1.5

Following [52], the material properties of the tissue model have been given in
Table I where the thickness of the corresponding layer is denoted by $d$ in
mm. The electrical conductivity, relative permittivity and mass density of the
layers are denoted by $\rho$, $\epsilon_{r}$ and $\sigma$ in Table I. Since
the tissue layer is placed with very close electrical distance
($0.016\lambda_{0}$, $\lambda_{0}$ being the free-space wavelength at $f=2.4$
GHz) from the antenna element, it lies in the radiative near field region at
$f=2.4$ GHz.

(a) Ref.,skin

(b) Ref.,fat

(c) X-slot,skin

(d) X-slot,fat

(e) Y-slot,skin

(f) Y-slot,fat

Figure 9: SAR distribution for three probe-fed antenna topologies with two
different tissue layers at $f=2.4$ GHz. Tissue is lying along the XZ-plane
with the dimension of $85$ mm $\times 60$ mm.

Being placed at horizontal separation, the tissue will be coupled by the
horizontal radiation components. The electromagnetic exposure of the tissue
layer can be characterised by specific absorption rate (SAR) parameter of
(10). The incident field on the tissue is directly related to the induced
characteristic modes of the antenna which vary its nature due to orientation
of the narrow slot. The field emerging out of the radiating patch antenna will
experience multiple reflection and scattering from the tissue layer. As a
consequence, the SAR distribution on the tissue will be perturbed accordingly.
Apart from the induced field $\vec{E}_{o}$, effective amount of the SAR value
is also dependent on $\rho$ and $\sigma$. For the case study, two different
tissues-skin and fat have been considered. The SAR distribution for the probe-
fed reference patch have been computed using Ansys HFSS and shown in Figs.
9(a)–9(f). It can be noted from the comparison of Figs. 9(a), 9(c), and 9(e)
that the X-directed slot reduces the electromagnetic energy absorption
compared to the reference antenna with no slot. The reason will be obvious
from the preceding current distribution and near filed distribution where the
Y-directed slot suppresses the Y-polarised current distribution of the
reference patch.

TABLE II: Maximum SAR value comparison of the probe-fed patch

Tissue | $\frac{max(SAR_{ref})}{max(SAR_{X-slot})}$ | $\frac{max(SAR_{ref})}{max(SAR_{Y-slot})}$
---|---|---
skin | $8.23$ | $0.8492$
fat | $7.77$ | $0.78$

Alternately, the Y-directed slot enhances the Y-polarised mode. As a
consequence, the SAR value is found to be slightly increased for the patch
with the Y-directed slot. The slot mode nature leads to generate two
concentric maxima regions in the skin layer of Fig. 9(c). Similar
characteristics can be found the examples with fat layer also. The SAR value
is reduced by the X-directed slot and enhanced by the Y-directed slot. Maximum
SAR values for the tissue layers have been compared in Table II for different
orientation of the slot. The next discussion considers microstrip line fed
patches for comparison.

Figure 10: Normalized distribution of the characteristic current modes at
$f=2.4$ GHz. (c) $1^{st}$ and (d) $2^{nd}$ modes of the plate with X-directed
slot.

###  III-B Microstrip line fed Patch

Similar patch elements have been considered but with microstrip line feeding
arrangement. Using Ansys Design Kit, the feed-line of the reference antenna is
tuned to 2.4 GHz with the same dimension of the radiating patch. The dimension
of the ground plane is considered to be 83.6 mm $\times$ 156 mm. The current
distribution have been compared in Figs. 10–10. Similar to the probe-fed
examples, the X-directed slot is found to suppress the Y-polarised edge mode
of Fig. 10. The impedance matching behaviour is also influenced by the
orientation of the slot in Fig. 10.

(a) $20log_{10}|\vec{E}_{1}|$ along X-axis

(b) $20log_{10}|\vec{E}_{1}|$ along Y-axis

Figure 11: Near field distribution due to the microstrip line-fed patch at
$f=2.4$ GHz.

(a) $\phi=0^{\circ}$-plane

(b) $\phi=0^{\circ}$-plane

Figure 12: Far field distribution due to the microstrip line-fed patch at
$f=2.4$ GHz.

(a) Ref.,skin

(b) Ref.,fat

(c) X-slot,skin

(d) X-slot,fat

(e) Y-slot,skin

(f) Y-slot,fat

Figure 13: SAR distribution for three line-fed antenna topologies with two
different tissue layers at $f=2.4$ GHz. Tissue is lying along the XZ-plane
with the dimension of $85$ mm $\times 40$ mm.

As shown in Figs. 11(a) and 11(b), the near field due to the horizontal
radiation is also decreased by the orthogonal slot. Henceforth, the intuitive
hints provided by the characteristic mode analysis is valid for both types of
feeding arrangements.

TABLE III: Maximum SAR value comparison of the line-fed patch

Tissue | $\frac{max(SAR_{ref})}{max(SAR_{X-slot})}$ | $\frac{max(SAR_{ref})}{max(SAR_{Y-slot})}$
---|---|---
skin | $418.3$ | $0.79$
fat | $336.56$ | $2.36$

The impacts of the current perturbation on the resultant SAR distribution have
been studied in Fig.13(a)–13(f) for both skin and fat tissue layers. As shown
in Fig.13(c) and 13(d), the X-directed slot reduces the SAR value compared to
the reference antennas of Fig.13(a) and 13(b). The nature of distribution is
also changed as there appears a concentric maxima around the centre for the
X-directed slot. However, for the Y-directed slot, other constitutive
parameters like $\sigma$ and $\rho$ also control the variation of SAR
following (10). As a consequence, the SAR value is enhanced by the Y-slot for
the skin layer where as for the fat layer, it is slightly increased.
Comparison of the maximum SAR values of the line fed antennas have been
reported in Table III. The impacts of X-directed slot is found to be maximum
for the microstrip line feeding with respect to the probe feeding because the
direction of the current in the feed line is parallel to the target
Y-polarised edge mode of the antenna. So the modal weighting coefficient
$\alpha_{n}$ of (3) starts to be an influencing parameter. As two vectors
$[V]$ and $[I_{m}]$ becomes parallel, their inner product holds higher value
in the line fed antennas. Thus, both the radiating element’s modal
characteristics and the orientation of the excitation signal control the
resultant radiation behaviour.

Determine the feed location from the maxima of the desired $i^{th}$ mode of
the reference patch. Determine the excitation vector $[V]$ based on the
corresponding feed location. Determine the modal coefficients $|\alpha_{n}|$
using the computed $[V]$. Compute the $1^{st}U$ number of significant modes
having significant interference to the desired $i^{th}$ mode based on the
higher value of $|\alpha_{n}|$, $n\neq i$. Determine the resultant current
density vector $\vec{J}_{und}=\sum_{n=1}^{U}|\alpha_{n}|\vec{J}_{n}$ where
$n\neq i$. Align the slot with the larger axis symmetrically orthogonal to the
direction of $\vec{J}_{und}$ in the centre of the patch. Figure 14: Procedure
for determining the orientation of the narrow slot.

In general, the characteristic mode analysis is found to provide systematic
guideline for maintaining the radiated signal’s polarization requirement and
SAR distribution along a desired direction. For arbitrary orientation of the
current distribution, the procedure of slot loading has been included in Fig.
14.

##  IV Conclusion

This work explored the impacts of slot loading in rectangular patch antennas.
Initially characteristic mode analysis of the radiating structures were
performed to understand the perturbation in the current and field
distribution. Later full-wave study had been carried out to investigate the
loading effects on the electromagnetic energy absorption by biological tissues
placed in the near field region. It was found that the prior knowledge of the
characteristic modes of the radiating geometry can provide suitable
information about the arrangement of the narrow slot for reducing the SAR
value of the neighbouring tissues. The orientation of the slot is found to
determine the amount of maximum radiation in a particular direction. Such
important feature of the CMA can be further utilised for focussing the power
of the wireless power transfer system to a target region of interest.

## References

  * [1] H. Elayan, R. M. Shubair, J. M. Jornet, A. Kiourti, and R. Mittra, “Graphene-Based Spiral Nanoantenna for Intrabody Communication at Terahertz,” in _2018 IEEE International Symposium on Antennas and Propagation USNC/URSI National Radio Science Meeting_ , Jul. 2018, pp. 799–800.
  * [2] H. Elayan, R. M. Shubair, and A. Kiourti, “On graphene-based THz plasmonic nano-antennas,” in _2016 16th Mediterranean Microwave Symposium (MMS)_ , Nov. 2016, pp. 1–3.
  * [3] R. M. Shubair and Y. L. Chow, “A simple and accurate approach to model the coupling of vertical and horizontal dipoles in layered media,” in _IEEE Antennas and Propagation Society International Symposium 1992 Digest_ , Jun. 1992, pp. 2309–2312 vol.4.
  * [4] A. A. Ibrahim and R. M. Shubair, “Reconfigurable band-notched UWB antenna for cognitive radio applications,” in _2016 16th Mediterranean Microwave Symposium (MMS)_ , Nov. 2016, pp. 1–4.
  * [5] M. S. Khan, A. Capobianco, S. M. Asif, A. Iftikhar, B. D. Braaten, and R. M. Shubair, “A pattern reconfigurable printed patch antenna,” in _2016 IEEE International Symposium on Antennas and Propagation (APSURSI)_ , Jun. 2016, pp. 2149–2150.
  * [6] R. M. Shubair, A. Salah, and A. K. Abbas, “Novel implantable miniaturized circular microstrip antenna for biomedical telemetry,” in _2015 IEEE International Symposium on Antennas and Propagation USNC/URSI National Radio Science Meeting_ , Jul. 2015, pp. 947–948.
  * [7] S. Alharbi, R. M. Shubair, and A. Kiourti, “Flexible Antennas for Wearable Applications: Recent Advances and Design Challenges,” in _12th European Conference on Antennas and Propagation (EuCAP 2018)_ , Jan. 2018.
  * [8] O. M. Khan, Q. U. Islam, R. M. Shubair, and A. Kiourti, “Novel multiband Flamenco fractal antenna for wearable WBAN off-body communication applications,” in _2018 International Applied Computational Electromagnetics Society Symposium (ACES)_ , Mar. 2018, pp. 1–2.
  * [9] M. S. Khan, A. Capobianco, S. M. Asif, A. Iftikhar, B. D. Braaten, and R. M. Shubair, “A properties comparison between copper and graphene-based UWB MIMO planar antennas,” in _2016 IEEE International Symposium on Antennas and Propagation (APSURSI)_ , Jun. 2016, pp. 1767–1768.
  * [10] S. Ghosal, A. De, A. Chakrabarty, and R. M. Shubair, “Characteristic Mode Analysis of Slot Loading in Microstrip Patch Antenna,” in _2018 IEEE International Symposium on Antennas and Propagation USNC/URSI National Radio Science Meeting_ , Jul. 2018, pp. 1523–1524.
  * [11] O. M. Khan, R. M. Shubair, and Q. U. Islam, “Second order flamenco fractal antenna for industrial scientific and medical applications,” in _2017 International Conference on Electrical and Computing Technologies and Applications (ICECTA)_ , Nov. 2017, pp. 1–3.
  * [12] M. S. Khan, A. Iftikhar, A.-D. Capobianco, R. M. Shubair, and B. Ijaz, “Pattern and frequency reconfiguration of patch antenna using PIN diodes,” _Microwave and Optical Technology Letters_ , vol. 59, no. 9, pp. 2180–2185, 2017. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/mop.30709
  * [13] R. M. Shubair and Y. L. Chow, “Complex images of vertical antenna penetrating a dielectric half-space,” in _Proceedings of IEEE Antennas and Propagation Society International Symposium_ , Jun. 1993, pp. 1187–1190 vol.3.
  * [14] R. A. Kipp and C. H. Chan, “A numerically efficient technique for the method of moments solution for planar periodic structures in layered media,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 42, no. 4, pp. 635–643, Apr. 1994.
  * [15] O. M. Khan, R. M. Shubair, Q. U. Islam, and I. Rashid, “Triband Metamaterial Embedded Implantable Antenna for Biotelemetry Applications,” in _2018 IEEE International Symposium on Antennas and Propagation USNC/URSI National Radio Science Meeting_ , Jul. 2018, pp. 213–214.
  * [16] R. Karli, H. Ammor, R. M. Shubair, M. I. A. Hajri, and A. Hakam, “Minituarized printed L-shaped antenna with parasitic element for ultra-wideband applicatons,” in _2016 17th International Symposium on Antenna Technology and Applied Electromagnetics (ANTEM)_ , Jul. 2016, pp. 1–2.
  * [17] R. M. Shubair, “An optimized complex-image prediction model for a closed-form solution of vertical antennas above a dielectric half-space,” in _1998 Symposium on Antenna Technology and Applied Electromagnetics_ , Aug. 1998, pp. 283–286.
  * [18] R. M. Shubair and Y. L. Chow, “A technique for the efficient computation of the periodic Green’s function in layered dielectric media,” in _1992 IEEE MTT-S Microwave Symposium Digest_ , Jun. 1992, pp. 393–396 vol.1.
  * [19] R. M. Shubair, Y. L. Chow, and J. J. Yang, “A full wave complex image solution of horizontal electric dipoles in a thick dielectric sheet,” in _Antennas and Propagation Society Symposium 1991 Digest_ , Jun. 1991, pp. 1372–1375 vol.3.
  * [20] R. M. Shubair and Y. L. Chow, “Analysis of the coupling of vertical and horizontal dipoles in layered media,” in _1992 Symposium on Antenna Technology and Applied Electromagnetics_ , Aug. 1992, pp. 211–216.
  * [21] ——, “A combined complex image-induced EMF solution of vertical antennas above a dielectric half-space,” in _Proceedings of IEEE Antennas and Propagation Society International Symposium_ , Jun. 1993, pp. 1030–1033 vol.2.
  * [22] ——, “Complex images for the calculation of the electrostatic potential in multilayered media,” in _Proceedings of IEEE Antennas and Propagation Society International Symposium_ , Jun. 1993, pp. 166–169 vol.1.
  * [23] R. Karli, H. Ammor, R. M. Shubair, M. I. AlHajri, R. Alkurd, and A. Hakam, “Miniature planar ultra-wide-band microstrip antenna for breast cancer detection,” in _2016 16th Mediterranean Microwave Symposium (MMS)_ , Nov. 2016, pp. 1–4.
  * [24] R. M. Shubair and Y. Chow, “A simple and accurate complex image interpretation of vertical antennas present in contiguous dielectric half-spaces,” _IEEE transactions on antennas and propagation_ , vol. 41, no. 6, pp. 806–812, 1993.
  * [25] ——, “Efficient computation of the periodic green’s function in layered dielectric media,” _IEEE transactions on microwave theory and techniques_ , vol. 41, no. 3, pp. 498–502, 1993.
  * [26] M. S. Khan, A.-D. Capobianco, S. M. Asif, D. E. Anagnostou, R. M. Shubair, and B. D. Braaten, “A compact csrr-enabled uwb diversity antenna,” _IEEE Antennas and Wireless Propagation Letters_ , vol. 16, pp. 808–812, 2016.
  * [27] R. Shubair and Y. Chow, “A closed-form solution of vertical dipole antennas above a dielectric half-space,” _IEEE transactions on antennas and propagation_ , vol. 41, no. 12, pp. 1737–1741, 1993.
  * [28] M. Al-Nuaimi, R. Shubair, and K. Al-Midfa, “Direction of arrival estimation in wireless mobile communications using minimum variance distortionless response,” in _The Second International Conference on Innovations in Information Technology (IIT’05)_ , 2005, pp. 1–5.
  * [29] G. Nwalozie, V. Okorogu, S. Maduadichie, and A. Adenola, “A simple comparative evaluation of adaptive beam forming algorithms,” _International Journal of Engineering and Innovative Technology (IJEIT)_ , vol. 2, no. 7, 2013.
  * [30] W. Che, C. Li, P. Russer, and Y. Chow, “Propagation and band broadening effect of planar integrated ridged waveguide in multilayer dielectric substrates,” in _2008 IEEE MTT-S International Microwave Symposium Digest_. IEEE, 2008, pp. 217–220.
  * [31] R. Shubair, “Robust adaptive beamforming using lms algorithm with smi initialization,” in _2005 IEEE Antennas and Propagation Society International Symposium_ , vol. 4. IEEE, 2005, pp. 2–5.
  * [32] A. Omar and R. Shubair, “Uwb coplanar waveguide-fed-coplanar strips spiral antenna,” in _2016 10th European Conference on Antennas and Propagation (EuCAP)_. IEEE, 2016, pp. 1–2.
  * [33] M. Bakhar and D. P. Hunagund, “Eigen structure based direction of arrival estimation algorithms for smart antenna systems,” _IJCSNS International Journal of Computer Science and Network Security_ , vol. 9, no. 11, pp. 96–100, 2009.
  * [34] M. S. Khan, A.-D. Capobianco, A. Iftikhar, R. M. Shubair, D. E. Anagnostou, and B. D. Braaten, “Ultra-compact dual-polarised uwb mimo antenna with meandered feeding lines,” _IET Microwaves, Antennas & Propagation_, vol. 11, no. 7, pp. 997–1002, 2017.
  * [35] R. M. Shubair, A. M. AlShamsi, K. Khalaf, and A. Kiourti, “Novel miniature wearable microstrip antennas for ism-band biomedical telemetry,” in _2015 Loughborough Antennas & Propagation Conference (LAPC)_. IEEE, 2015, pp. 1–4.
  * [36] M. Y. ElSalamouny and R. M. Shubair, “Novel design of compact low-profile multi-band microstrip antennas for medical applications,” in _2015 Loughborough Antennas & Propagation Conference (LAPC)_. IEEE, 2015, pp. 1–4.
  * [37] A. A. Ibrahim, J. Machac, and R. M. Shubair, “Compact uwb mimo antenna with pattern diversity and band rejection characteristics,” _Microwave and Optical Technology Letters_ , vol. 59, no. 6, pp. 1460–1464, 2017.
  * [38] R. Shubair and A. Merri, “Convergence of adaptive beamforming algorithms for wireless communications,” in _Proc. IEEE and IFIP International Conference on Wireless and Optical Communications Networks_ , 2005, pp. 6–8.
  * [39] J. An, G.-M. Wang, C.-X. Zhang, and H.-Y. Zeng, “Compact circularly polarized omnidirectional microstrip antenna,” _Microwave and Optical Technology Letters_ , vol. 51, no. 11, pp. 2643–2646, 2009.
  * [40] K. Ray, R. Kulkarni, and B. K. Ramkrishnan, “Doa estimation in a multipath environment using covariance differencing and iterative forward and backward spatial smoothing,” in _2008 International Conference on Recent Advances in Microwave Theory and Applications_. IEEE, 2008, pp. 794–796.
  * [41] A. Hakam, M. Hussein, M. Ouda, R. Shubair, and E. Serria, “Novel circular antenna with elliptical rings for ultra-wide-band,” in _2016 10th European Conference on Antennas and Propagation (EuCAP)_. IEEE, 2016, pp. 1–4.
  * [42] M. I. Hussein, A. Hakam, M. Ouda, and R. M. Shubair, “Compact low-profile planar elliptical antenna for uwb applications,” in _2016 10th European Conference on Antennas and Propagation (EuCAP)_. IEEE, 2016, pp. 1–2.
  * [43] E. M. Al-Ardi, R. M. Shubair, and M. E. Al-Mualla, “Investigation of high-resolution DOA estimation algorithms for optimal performance of smart antenna systems,” in _4th International Conference on 3G Mobile Communication Technologies_ , Jan. 2003, pp. 460–464.
  * [44] ——, “Performance evaluation of direction finding algorithms for adapative antenna arrays,” in _10th IEEE International Conference on Electronics, Circuits and Systems, 2003. ICECS 2003. Proceedings of the 2003_ , vol. 2, Dec. 2003, pp. 735–738 Vol.2.
  * [45] J. Wiart, C. Dale, A. V. Bosisio, and A. Le Cornec, “Analysis of the influence of the power control and discontinuous transmission on rf exposure with gsm mobile phones,” _IEEE Transactions on electromagnetic compatibility_ , vol. 42, no. 4, pp. 376–385, 2000.
  * [46] P. Bernardi, M. Cavagnaro, S. Pisa, and E. Piuzzi, “Power absorption and temperature elevations induced in the human head by a dual-band monopole-helix antenna phone,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 49, no. 12, pp. 2539–2546, 2001.
  * [47] A. Gati, A. Hadjem, M.-F. Wong, and J. Wiart, “Exposure induced by wcdma mobiles phones in operating networks,” _IEEE Transactions on wireless communications_ , vol. 8, no. 12, pp. 5723–5727, 2009.
  * [48] S. Kühn and N. Kuster, “Field evaluation of the human exposure from multiband, multisystem mobile phones,” _IEEE Transactions on Electromagnetic Compatibility_ , vol. 55, no. 2, pp. 275–287, 2012.
  * [49] B. Xu, M. Gustafsson, S. Shi, K. Zhao, Z. Ying, and S. He, “Radio frequency exposure compliance of multiple antennas for cellular equipment based on semidefinite relaxation,” _IEEE Transactions on Electromagnetic Compatibility_ , vol. 61, no. 2, pp. 327–336, 2018.
  * [50] D. G. Lopez, M. Ignatenko, and D. S. Filipovic, “Eigenmode prediction of high rf exposure frequency region inside vehicles,” _IEEE Transactions on Electromagnetic Compatibility_ , vol. 59, no. 1, pp. 43–47, 2016.
  * [51] S. Ghosal, A. De, A. P. Duffy, and A. Chakrabarty, “Selection of dominant characteristic modes,” _IEEE Trans Electromagn Compat_ , pp. 1–10, 2019.
  * [52] X.-Q. Zhu, Y.-X. Guo, and W. Wu, “Miniaturized dual-band and dual-polarized antenna for mban applications,” _IEEE Trans. Antennas Propagat._ , vol. 64, no. 7, pp. 2805–2814, 2016.
  * [53] B. J. DeLong, A. Kiourti, and J. L. Volakis, “A radiating near-field patch rectenna for wireless power transfer to medical implants at 2.4 ghz,” _IEEE Journal of Electromagnetics, RF and Microwaves in Medicine and Biology_ , vol. 2, no. 1, pp. 64–69, 2018.
  * [54] A. S. Poon, S. O’Driscoll, and T. H. Meng, “Optimal frequency for wireless power transmission into dispersive tissue,” _IEEE Trans. Antennas Propagat._ , vol. 58, no. 5, pp. 1739–1750, 2010.
  * [55] A. S. Poon, “Electromagnetic field focusing for short-range wireless power transmission,” in _Radio and Wireless Symposium (RWS), 2012 IEEE_. IEEE, 2012, pp. 115–118.
  * [56] G. G. Bellizzi, L. Crocco, G. M. Battaglia, and T. Isernia, “Multi-frequency constrained sar focusing for patient specific hyperthermia treatment,” _IEEE Journal of Electromagnetics, RF and Microwaves in Medicine and Biology_ , vol. 1, no. 2, pp. 74–80, 2017.
  * [57] Y. Fan, J. Huang, T. Chang, and X. Liu, “A miniaturized four-element mimo antenna with ebg for implantable medical devices,” _IEEE Journal of Electromagnetics, RF and Microwaves in Medicine and Biology_ , 2018.
  * [58] S. Ghosal and S. R. B. Chaudhuri, “Analysis of a rectangular slot on a microstrip patch antenna with an equivalent circuit model,” in _Applied Electromagnetics Conference (AEMC), 2013 IEEE_. IEEE, 2013, pp. 1–2.
  * [59] M. U. Khan, M. S. Sharawi, and R. Mittra, “Microstrip patch antenna miniaturisation techniques: a review,” _IET Microw., Antennas Propag._ , vol. 9, no. 9, pp. 913–922, 2015.
  * [60] J.-Y. Sze and K.-L. Wong, “Slotted rectangular microstrip antenna for bandwidth enhancement,” _IEEE Trans. Antennas Propagat._ , vol. 48, no. 8, pp. 1149–1152, 2000.
  * [61] P. Juyal and L. Shafai, “Sidelobe reduction of tm 12 mode of circular patch via nonresonant narrow slot,” _IEEE Trans. Antennas Propagat._ , vol. 64, no. 8, pp. 3361–3369, 2016.
  * [62] D. Pozar, “A reciprocity method of analysis for printed slot and slot-coupled microstrip antennas,” _IEEE Trans. Antennas Propagat._ , vol. 34, no. 12, pp. 1439–1446, 1986.
  * [63] S. Ghosal, A. De, A. Chakrabarty, and R. M. Shubair, “Characteristic mode analysis of slot loading in microstrip patch antenna,” _IEEE APS/URSI Symp. 2018, Boston, USA_ , 2018.
  * [64] ——, “Analysis of slot loading in elliptical patch- a characteristic mode approach,” _IEEE ANTEM 2018, Canada_ , 2018.
  * [65] R. Harrington and J. Mautz, “Theory of characteristic modes for conducting bodies,” _IEEE Trans. Antennas Propag_ , vol. 19, no. 5, pp. 622–628, 1971\.
  * [66] F. H. Lin and Z. N. Chen, “A method of suppressing higher-order modes for improving radiation performance of metasurface multiport antennas using characteristic mode analysis,” _IEEE Trans. Antennas Propagat._ , 2018.
  * [67] J. Zhao, Y. Chen, and S. Yang, “In-band radar cross section reduction of slot antenna using characteristic modes,” _IEEE Antennas Wireless Propag. Lett._ , 2018.
  * [68] U. G. Matlab, “The mathworks,” _Inc., Natick, MA_ , vol. 1992, 1760.
  * [69] S. N. Makarov, _Antenna and EM Modeling with MATLAB_. Wiley-Interscience,, 2002.
  * [70] H. ANSYS, “Ansoft corp., pittsburgh, pa,” 2015.


# Harmonization of Medical Data with Flow-based Causal Inference

Rongguang Wang  Pratik Chaudhari  Christos Davatzikos

###### Abstract

Heterogeneity in medical data, e.g., caused due to the diversity of data
collected at different sites of a clinical study, is a fundamental hurdle for
accuracy and generalization of machine learning predictive models. This paper
presents a normalizing-flow-based method to perform counterfactual inference
upon a structural causal model (SCM) to harmonize such data. We formulate a
causal model for observed effects (e.g., brain magnetic resonance imagining
data) that result from known confounders (e.g., site, gender and age) and
exogenous noise variables. Our method exploits the bijection induced by flow
for harmonization. We can infer the posterior of exogenous variables,
intervene on observations, and draw samples from the resultant SCM to obtain
the counterfactuals. We evaluate on multiple, large, real-world medical
datasets and existing causal inference datasets to observe that this method
leads to better cross-domain generalization compared to state-of-the-art
algorithms. Further experiments that evaluate the quality of confounder-
independent data generated by our model using regression and classification
tasks are provided.

##  1 Introduction

Start with example of multi-site data. Draw a plot of the comparison of
invariant risk minimization and fine-tuning. The former is insensitive to the
confounder but loses out on accuracy. The latter improves the accuracy but has
to use information from the confounding variable.

This paper focuses on harmonization of medical data. We model harmonization a
counterfactual query on a structural causal model. Through the paper, we
consider the case when the SCM is known, although there may exist exogenous
noise variables that model the effect of unknown causes.

The recent advances of deep learning also show its expressive capacity
function fitting ability in many applications such as computer-aided diagnosis
(esteva2017dermatologist) and segmentation (menze2014multimodal) in medical
imaging. However, these sufficiently trained models are sensitive to data
covariates changes. As shown in Figure 1 (a), the common cause $z$ of
observation $x$ and label $y$ is the confounding or confounder in statistics.
Usually, the likelihood $P(y|x)$ estimated by learning models ignores
confounding $z$, and this could lead to large generalization gap when testing
the model on data with confounding variations. In magnetic resonance imaging
(MRI), variations in scanner characteristics such as the magnetic field
strength, scanner vendor and acquisition protocols may have strong effects on
collected data. Therefore, it’s essential to tackle confoundings such as site
variations in machine learning to perform reliable analysis of large
populations across sites (davatzikos2019machine).

In this work, we tackle the counfounding problem from a causal perspective. By
explicitly modeling causal relationship of observations and the confoundings,
and acknowledging the difference between causation and correlation, we can
intervene the site variable and generate site-independent data. We propose a
structural causal equation to model the causal relationship between images and
its covariates such as sex, age, and site, which is based on a tractable
generative network, namely, normalizing flows (papamakarios2019normalizing).
This method has several advantages over existing methods. Firstly, the flow-
based structural causal model takes both known and unknown confoundings into
consideration, which preserves the nature of the original observation.
Secondly, since the transformations in flow models are invertible, they
provide exact posterior, whereas amortized models such as variational auto-
encoders (VAEs) (kingma2013auto) and generative adversarial networks (GANs)
(goodfellow2014generative) can only approximate the solution. Finally, the
flow models are more expressive in density estimation compared to VAEs and
GANs.

We demonstrate our proposed method on several large-scale brain imaging
datasets, showing substantial improvement compared to a recently published
baseline method on both age regression and Alzheimer’s disease classification
tasks. In particular, the testing classification accuracies of our models are
close to the ones trained directly on the target site.

###  1.1 Contributions

A formulation of the harmonization problem for brain MRI data using the causal
inference framework

A normalizing flow-based method for performing counterfactual queries on the
SCM

state of the art results on real-world and challenging medical datasets.
Thorough ablation studies…

##  2 Related work

To remove the undesired counfoundings, especially sites or scanners, from
imaging data, a wide range of contributions has been made by employing the
recent advances in statistical and machine learning (johnson2007adjusting;
pomponio2020harmonization; moyer2018invariant; moyer2020scanner;
robinson2020image; bashyam2020medical). Several prior works approach the
problem by removing site-related location-scale effects (johnson2007adjusting;
pomponio2020harmonization), extracting scanner-invariant representations
(moyer2018invariant; moyer2020scanner), and learning mappings between site
(robinson2020image; bashyam2020medical). Based on parametric empirical bayes
(morris1983parametric), ComBat methods (johnson2007adjusting;
pomponio2020harmonization) produce site-removed image features by performing
location (mean) and scale (variance) adjustments to the data. The linear model
estimates the location and scale differences in images features cross-site,
while preserving other confounders such as sex and age. Since ComBat
(johnson2007adjusting; pomponio2020harmonization) only considered covariates
sex and age in modelling, other unknown variations such as race and disease
are removed together with the site variable, which might leads to
unsatisfactory performance in downstream tasks, e.g. disease diagnosis. On the
other hand, generative deep learning models such as VAEs (kingma2013auto) and
GANs (goodfellow2014generative) have been used in many works
(moyer2018invariant; moyer2020scanner; robinson2020image; bashyam2020medical).
In order to disentangle the scanner-specific information from the images,
(moyer2018invariant; moyer2020scanner) propose to minimize the mutual
information between the site variable and image embedding in the latent space
of a VAE with conditional decoder. Then, the site-invariant representations
can be composed with a reference site variable to reconstruct images at the
same site. To avoid the blurriness problem of VAEs, unsupervised image-to-
image translation has been proposed to map scans either between two sites
(robinson2020image) or to a reference domain (bashyam2020medical) using
CycleGAN (zhu2017unpaired). Although these methods perform well empirically,
the GANs suffer from convergence issues, especially for 3D images.

##  3 Background

###  3.1 Ladder of Causation

(a) Confounding

(b) SCM

Figure 1: (a) Causal structure of $z$ confounding the effect of $x$ on $y$.
$x$ and $y$ are observations and labels separately. $z$ is the common cause of
$x$ and $y$, which is known as confounder or confounding; (b) Structural
causal model. $x_{k}$ is the observation, $\epsilon_{k}$ is the exogeneous
noise of $x_{k}$, and $\boldsymbol{pa}_{k}$ are the direct causes of $x_{k}$.
$f_{k}$ is an invertible transformation indicated by the bidirectional arrows.
The black and white arrowheads represent generative and abduction directions
separately. The edge ending with black dot denotes that $x_{k}$ is conditioned
on its direct causes $\boldsymbol{pa}_{k}$.

The causal mechanism proposed by Judea Pearl is described as the three-level
of causation ladder (pearl2009causality) which consists of association,
intervention, and counterfactuals. Association describes the correlations in
the data which reasoning about the observed data, which is the lower rung of
the ladder. This layer is characterized by conditional probability sentence
$P(y|x)$ which states that the probability of event $Y=y$ given we observed
event $X=x$. Intervention, at a higher level than association, involves not
just seeing what is, but changing what we see. It requires structural
assumptions based on prior knowledge about the underlying data generation. A
typical question at this level would be: _what happens if we do … ?_. We
usually find sentences of the type $P(y|do(\tilde{x}),z)$ denoting the
probability of event $Y=y$ given that we intervene and set the value of $X$ to
$\tilde{x}$ and subsequently observe confounder event $Z=z$. The highest rung
on the ladder of causation is counterfactuals. Counterfactuals provide
alternative outcomes from the generative process with retrospective queries. A
typical question in the counterfactual category is _what if we were to do
differently?_. This layer is characterized by expressions of the type
$P(\tilde{y}_{\tilde{x}}|x,y)$ which stand for the probability of event
$Y=\tilde{y}$ had $X$ been $\tilde{x}$, given that we actually observed $X$ to
be $x$ and and $Y$ to be $y$. Counterfactuals are placed at the top of the
hierarchy because they subsume interventional and associational questions. If
we have a model that can answer counterfactual queries, we can also answer
questions about interventions and observations. Counterfactuals expressions
can be possessed by structural causal models (SCMs) (peters2017elements).

###  3.2 Structural Causal Models (SCMs)

Consider a collection of $K$ random variables
$\boldsymbol{x}=(x_{1},\cdots,x_{K})$, a structural causal model
(peters2017elements; scholkopf2019causality)
$\mathcal{M}:=(\boldsymbol{S},P(\boldsymbol{\epsilon}))$ consists of a
collection $\boldsymbol{S}=(f_{1},\cdots,f_{K})$ of structural assignments
$x_{k}:=f_{k}(\epsilon_{k};\boldsymbol{pa}_{k})$, where $\boldsymbol{pa}_{k}$
is the set of parents (direct causes) of $x_{k}$, and a joint distribution
$P(\boldsymbol{\epsilon})=\prod^{K}_{k=1}P(\epsilon_{k})$ over mutually
independent exogenous noises $\boldsymbol{\epsilon}$. The exogenous noises
$\bm{\epsilon}$ are associated with the unknown sources of variation. The
model can be represented by a directed acyclic graph (DAG) with edges pointing
from causes to effects, which is a probabilistic graphical model (PGM) with
causal relations as shown in Figure 1 (b). Since each variable $\bm{x}$ is
independent of its non-effects given its direct causes (also known as causal
Markov condition), the joint observational distribution of SCM can be
factorized as

$P_{\mathcal{M}}(\bm{x})=\prod^{K}_{k=1}P_{\mathcal{M}}(x_{k}|\bm{pa}_{k}),$
(1)

where each conditional distribution $P_{\mathcal{M}}(x_{k}|\bm{pa}_{k})$ is
determined by the corresponding structural assignment and noise distribution
(pearl2009causality).

###  3.3 Counterfactual Inference

Given SCMs, the counterfactual query is formulated as a three-step process,
namely, abduction, action, and prediction (pearl2009causality;
peters2017elements; pearl2009causal). Firstly, we predict the exogenous noise
$\bm{\epsilon}$ based on the observations $\bm{x}$, which means we inferring
the posterior noise distribution:

$P_{\mathcal{M}}(\boldsymbol{\epsilon}|\boldsymbol{x})=\prod^{K}_{k=1}P_{\mathcal{M}}(\epsilon_{k}|x_{k},\boldsymbol{pa}_{k}).$
(2)

Then, we perform an intervention by replacing the structural assignment(s)
denoted as $do(x_{k}:=\tilde{x}_{k})$ which makes the variable independent of
its former parents $\boldsymbol{pa}_{k}$. Note that this will result in a
modified SCM:

$\mathcal{\tilde{M}}=\mathcal{M}_{\boldsymbol{x};do(\tilde{x}_{k})}=(\tilde{\boldsymbol{S}},P_{\mathcal{M}}(\boldsymbol{\epsilon}|\boldsymbol{x})),$
(3)

where each intervened variable $x_{k}$ is independent of its former parents
$\boldsymbol{pa}_{k}$ and of its exogenous noise $\epsilon_{k}$. Now, the
noise distribution in the modified SCM is the posterior
$P_{\mathcal{M}}(\boldsymbol{\epsilon}|\boldsymbol{x})$ obtained in the
abduction step. Finally, we predict the counterfactuals by sampling from the
distribution $P_{\mathcal{\tilde{M}}}(\boldsymbol{x})$ entailed by the
modified SCM.

###  3.4 Normalizing flows

To utilize the powerful function approximation ability of deep learning, we
model the structural assignments with neural networks. However, the abduction
step in counterfactual inference is potentially intractable due to posterior
noise distribution $P_{\mathcal{M}}(\boldsymbol{\epsilon}|\boldsymbol{x})$
computation. To enable tractable counterfactual inference, we use normalizing
flows to model the mechanism. Normalizing flows are used as density estimator
which modeling complex probability distributions using a set of
transformations from simple base distributions in the same dimensionality
(papamakarios2019normalizing; papamakarios2017masked;
dolatabadi2020invertible; durkan2019neural).

Given observed variables $\boldsymbol{x}$, and base variables
$\boldsymbol{\epsilon}\sim p(\boldsymbol{\epsilon})$, we would like to find
diffeomorphism, i.e., an invertible and differentiable, transformation $f$
such that $\boldsymbol{x}=f(\boldsymbol{\epsilon})$. The output distribution
$p(\boldsymbol{x})$ can be computed by a change of variables:

$p(\boldsymbol{x})=p(\boldsymbol{\epsilon})|\det\nabla
f(\boldsymbol{\epsilon})|^{-1},\
\boldsymbol{\epsilon}=f^{-1}(\boldsymbol{x}).$ (4)

Specifically, the function $f$ compresses and expands the density of base
distribution $p(\boldsymbol{\epsilon})$, which is quantified by the
determinant of the Jacobian of the transformation. We usually chose a simple
base distribution, such as a standard normal. Since diffeomorphic
transformations are composable, we can build complex transformations by
composing multiple instances of simple transformations, without compromising
the requirement of invertibility and differentiability. Hence, we can still
calculate the density $p(\boldsymbol{x})$ using Equation 4. The transformation
$f$ and the base distribution $p(\boldsymbol{\epsilon})$ can be parameterized
with neural networks. Given a dataset
$\mathcal{D}=\\{\boldsymbol{x}^{i}\\}^{N}_{i=1}$, the flow model can be
optimized by maximizing the total log-likelihood $\sum_{i}\log
p(\boldsymbol{x}^{i})$ with respect to the parameters of transformation $f$.

Conditional normalizing flows. Extending normalizing flows to conditional
densities, we parameterize the transformation as
$x=f(\epsilon;\boldsymbol{pa}_{X})$. Then, the output density of the
conditional flows is

$p(x|\boldsymbol{pa}_{X})=p(\epsilon)\cdot|\det\nabla_{\epsilon}f(\epsilon;\boldsymbol{pa}_{X})|^{-1},\
\epsilon=f^{-1}(x;\boldsymbol{pa}_{X}).$ (5)

Benefiting from the invertible mechanism, the noise posterior can be obtained
deterministically and exactly by evaluating at
$\epsilon=f^{-1}(x;\boldsymbol{pa}_{X})$.

##  4 Method

###  4.1 Algorithm

1: Initialize: flow parameters $\theta$ and learning rate $\alpha$

2: for each training iteration $i$ do

3: $(x^{i},\boldsymbol{pa}_{X}^{i})\in(x,\boldsymbol{pa}_{X})$

4: $\mathcal{L}(\theta)=\log
p(f_{\theta}^{-1}(x^{i};\boldsymbol{pa}_{X}^{i}))-\log|\det
J_{f^{-1}_{\theta}}(x^{i};\boldsymbol{pa}_{X}^{i})|$

5: $\theta\leftarrow\theta-\alpha\nabla_{\theta}\mathcal{L}(\theta)$

6: end for

7: Inference:
$P_{\mathcal{M}}(\boldsymbol{\epsilon}|\boldsymbol{x})=\prod^{K}_{k=1}p(f_{\theta}^{-1}(x_{k}),\boldsymbol{pa}_{k})$

8: Action $do(\cdot)$: $x_{k}:=\tilde{x}_{k}$ or
$x_{k}:=\tilde{f}_{k}(\epsilon_{k};\tilde{\boldsymbol{pa}}_{k})$

9: Sampling: $\boldsymbol{x}^{h}\sim P_{\mathcal{\tilde{M}}}(\boldsymbol{x})$
where
$\mathcal{\tilde{M}}=(\tilde{\boldsymbol{S}},P_{\mathcal{M}}(\boldsymbol{\epsilon}|\boldsymbol{x}))$

10: Return: $\boldsymbol{x}^{h}$

Algorithm 1 Harmonization with flow-based SCM.

###  4.2 SCM for brain MRI data

Explain the SCM

Explain how the conditional flow is used to perform the action step

Figure 2: Computational graph of the structural causal model for brain
imaging. We denote brain ROI volumes as $x$, sex as $s$, age as $a$, and
imaging site as $t$. The white dot indicates the information aggregation of
the direct causes $s$, $a$, and $t$.

In this section, we firstly train the flow-based SCM with passively observed
data (image features $x$, sex $s$, age $a$, and site $t$) from iSTAGING and
ADNI datasets, where we fitting the flow model with observations. Then, we
replace the structural assignment for site $t$ variable to a fixed one to make
sure that all image features $x$ are conditioned on the same site $t$. This is
the intervention step which builds on the structural assumption that sex $s$,
age $a$, and site $t$ are the causes of the image features $x$. Finally, we
sample from the modified flow-based SCM to obtain site-independent image
features. This step is achieved by performing counterfactural inference where
posterior distributions of exogeneous noises are computed by applying the
inverse transformation of the flow model to the observations.

We construct the flow-based SCM with image features (ROI volumes), clinical
variables, such as sex, age, and scanner information site where sex, age and
site are the direct causes (parents) of the brain image features as shown in
Figure 2. In the experiments, we assume the parent variables ($s$, $a$, and
$t$) are mutually independent and the image features $x$ are conditioned on
the aggregation of their direct causes.

###  4.3 An example with synthetic data

This is the example that you had created with a known causal model?

##  5 Experiments

In this section, we firstly introduce several large-scale real-world brain
imaging datasets for experiments. Then, we consider the problem of removing
confounding (scanner effects) in the brain images cross-site by modelling the
structural causal model and performing counterfactual inference. We evaluate
the performance of the propose method on two tasks: age regression and
Alzheimer’s disease classification. Finally, we show implementation details of
our flow-based model and the models used for regression and classification
tasks.

###  5.1 Datasets and Setup

Table 1: Summary of participant demographics in iSTAGING consortium. Age is described in format: mean $\pm$ std [min, max]. F and M in gender represent female and male separately. Field indicates the magnetic strength of the MRI scanners. Study | Subject | Age | Gender (F/M) | Field
---|---|---|---|---
BLSA-1.5T | 157 | 69.1 $\pm$ 8.5 [48.0, 85.0] | 66 / 91 | 1.5T
BLSA-3T | 960 | 65.0 $\pm$ 14.7 [22.0, 93.0] | 525 / 435 | 3T
UKBB | 2202 | 62.8 $\pm$ 7.3 [45.0, 79.0] | 1189 / 1013 | 3T
SHIP | 2739 | 52.6 $\pm$ 13.7 [21.2, 90.4] | 1491 / 1248 | 1.5T
Table 2: Summary of participant demographics in ADNI dataset. Age is described in format: mean $\pm$ std [min, max]. F and M in gender represent female and male separately. Field indicates the magnetic strength of the MRI scanners. Study | Subject | CN | AD | Age | Gender (F/M) | Field
---|---|---|---|---|---|---
ADNI-1 | 422 | 229 | 193 | 75.5 $\pm$ 6.2 [55.0, 90.9] | 201 / 221 | 1.5T
ADNI-2/GO | 441 | 294 | 147 | 73.4 $\pm$ 6.8 [55.4, 90.3] | 221 / 220 | 3T

We evaluate the proposed method on 3D T1-weighted brain magnetic resonance
imaging (MRI) scans from two consortia: Alzheimer’s Disease Neuroimaging
Initiative (ADNI) (jack2008alzheimer) study and iSTAGING consortium. We use
two phases (ADNI-1 and ADNI-2/GO) in ADNI dataset and four studies including
Baltimore Longitudinal Study of Aging (BLSA) (resnick2003longitudinal;
armstrong2019predictors) 1.5T SPGR and 3T MPRAGE studies, Study of Health in
Pomerania (SHIP) (hegenscheid2009whole), and UK Biobank (UKBB) (sudlow2015uk)
study in iSTAGING consortium. Both consortia consist of participants in
different genders and wide range of age. The brain scans are acquired from
multiple scanners with different scale of magnetic strength and acquisition
protocols. The detailed demographic information of iSTAGING and ADNI datasets
are shown in Table 1, and Table 2 separately. In the experiments, we perform
age prediction task using iSTAGING consortium where the age of participants is
range from $21$ to $93$ years old. We also demonstrate the effectiveness of
our method on Alzheimer’s disease (AD) classification task which is more
challenging than age prediction using ADNI dataset where the diagnosis groups
including cognitive normal (CN) and AD.

We perform a pipeline of preprocessing steps, including bias-filed correction
(tustison2010n4itk), brain tissue extraction via skull-stripping
(doshi2013multi), and multi-atlas segmentation (doshi2016muse), on the images.
In the last step, each scan is segmented into 145 anatomical regions of
interests (ROIs) spanning the entire brain. Then, we calculate the volume of
each ROI as image feature. The list of ROIs are provided in Appendix
LABEL:app:roi_list.

###  5.2 Details of training

We examine the three level of causation ladder hierarchy by removing the site
effects in iSTAGING and ADNI datasets. We use normalizing flows for the
structural assignments of all variables except sex and site. We directly learn
the binary probability of sex $s$ and categorical probability of site $t$. The
mechanism is modelled as follows:

$\begin{split}&s:=\epsilon_{S},\\\ &a:=f_{A}(\epsilon_{A}),\\\
&t:=\epsilon_{T},\\\
&x:=f_{X}(\epsilon_{X};s,a,t),\qquad\quad\end{split}\quad\begin{split}&\epsilon_{S}\sim\text{Ber}(p_{\theta}^{S});\\\
&\epsilon_{A}\sim\mathcal{N}(0,1);\\\
&\epsilon_{T}\sim\text{Cat}(K,p_{\theta}^{T});\\\
&\epsilon_{X}\sim\mathcal{N}(0,1),\\\ \end{split}$ (6)

where $p_{\theta}^{S}$ and $p_{\theta}^{T}$ are the learnable mass functions
of the categorical distribution for variables sex $s$ and site $t$, and $K$ is
the number of site $t$. The modules indicated with $\theta$ are parameterized
using neural networks. We constrain age $a$ variable with lower bound
(exponential transform) and rescale it with fixed affine transform for
normalization. We use zero-mean unit-variance Gaussians as the base
distributions for all exogeneous noise distribution $P(\boldsymbol{\epsilon})$
if possible. $\text{Spline}_{\theta}$ transformation refers to the linear
neural spline flows (dolatabadi2020invertible). The
$\text{ConditionalTransform}_{\theta}(\cdot)$ can be conditional Affine or
conditional spline transform, which reparameterizes the noise distribution
into another Gaussian distribution. We use linear (dolatabadi2020invertible)
and quadratic (durkan2019neural) autoregressive neural spline flows for the
conditional spline transform, which are more expressive compared to the Affine
flows. The transformation parameters of the
$\text{ConditionalTransform}_{\theta}(\cdot)$ are predicted by a context
neural network taking $\cdot$ as input. The context networks are implemented
as fully-connected networks for Affine and Spline flows separately (details in
Appendix LABEL:app:flow_arch).

We implement the transformations in flow-based SCM with PyTorch
(paszke2019pytorch) and Pyro (bingham2019pyro). Adam (kingma2014adam) with
initial learning rate $3\text{e}-4$ and weight decay $1\text{e}-4$ is used for
optimization. We utilize a step learning rate scheduler with decay milestones
at $50\%$ and $75\%$ of the total epochs in training. The batch size of the
training process is $64$. We use 32 Monte Carlo samples for counterfacturals
estimation. All models are trained for at most $100$ epochs and we report
results from the model with the best validation log-likelihood.

The model of age regression and AD classification tasks is implemented with a
three-layer multi-layer perceptron (MLP). The network architecture is provided
in Appendix LABEL:app:mlp_arch. We use stochastic gradient descent (SGD) for
optimization with batch size $64$ and a learning rate of $1\text{e}-1$
(momentum $0.9$). The training would stop in at most $60$ epochs for all
models in the experiments. Huber loss and Binary cross-entropy are the loss
functions for age prediction and AD classification tasks separately. We use
mean absolute error (MAE) and classification accuracy as the evaluation
metrics for regression and classification respectively.

###  5.3 Baselines

To evaluate the performance of our proposed method, we compare with the state-
of-the-art algorithm ComBat (johnson2007adjusting; pomponio2020harmonization),
which is designed to remove the site effects in brain imaging. We implement
the ComBat algorithm with the publicly available code and apply the method to
the datasets following the best practice. We show the effectiveness of the
methods on both age prediction and AD classification tasks. The site-
removed/independent image features are expected to show consistent performance
on these tasks inter-site.

Invariant Risk Minimization.

Fine-tuning.

ComBat.

###  5.4 Density estimation

Figure 3: Comparison of normalized feature (hippocampus volume) distributions cross-site in iSTAGING consortium before and after apply the ComBat methods (ComBat-Linear and ComBat-GAM) and the proposed methods (SCM-Affine, SCM-LSpline, and SCM-QSpline). Table 3: Comparison of associative abilities of different type of flows on iSTAGING consortium and ADNI dataset. Study | Model | Log-likelihood
---|---|---
Synthetic | Affine | -1.5547
iSTAGING | Affine | 1.8817
Linear Spline | 17.2204
Quadratic Spline | 17.2397
ADNI | Affine | 1.8963
Linear Spline | 15.2715
Quadratic Spline | 15.2055

We quantitatively compare the associative capabilities of Affine, linear
Spline, and quadratic Spline flow models by evaluating their total log-
likelihood as shown in Table 3. For both iSTAGING and ADNI dataset, the
performance improves consistently with the model’s expressive power. The
Spline autoregressive flow models show much better density estimation ability
compared to simple Affine flow model, while quadratic Spline model shows
slightly higher log-likelihood than the linear Spline model in iSTAGING
dataset.

After we obtain the trained flow-based SCM, intervention and counterfactual
inference are performed for site-independent image features $x$ in iSTAGING
dataset. We show the feature (hippocampus volume) distributions of raw data,
ComBat (johnson2007adjusting; pomponio2020harmonization) transformed data, and
flow-based SCM generated data respectively in Figure 3. We find that the
feature distributions are not consistent with different means and variances
inter-site in raw data. The ComBat transformed feature distributions (linear
and GAM) show relatively consistent means but inconsistent variances resulting
from site location and scale effects removal process in the method. The flow-
based SCM generated data are site-independent where site effects are not
removed. Thus, we expect the feature distributions of the proposed method to
be closed to the raw data. As shown in the lower level of Figure 3, the Affine
model couldn’t capture the densities of the data, which agrees with the
unsatisfactory log-likelihood in Table 3. In contrast, both Spline (linear and
quadratic) models show matched feature distributions cross-site compared to
the raw data.

###  5.5 Age prediction

Table 4: Age prediction MAE comparison in iSTAGING consortium. Baseline show results trained on target site directly. The raw, ComBat, and counterfactuals columns use the raw features, features generated by the ComBat algorithms, and features generated by the proposed method. | Study | TarOnly | SrcOnly | Transfer | IRM | ComBat | Flow-based SCM (ours)
---|---|---|---|---|---|---|---
Linear | GAM | Affine | L-Spline | Q-Spline
Source | BLSA-3T | - | 10.4871 | 11.4589 | 11.4601 | 11.4696 | 11.4734 | 11.4713 | 10.6533 | 11.4689
Target | BLSA-1.5T | 6.7767 | 9.7854 | 6.4123 | 6.4898 | 6.6869 | 6.6203 | 6.6459 | 5.9314 | 6.7256
Target | UKBB | 6.1354 | 13.2142 | 6.3340 | 7.5689 | 6.5297 | 6.5883 | 6.5658 | 6.5413 | 6.4992
Target | SHIP | 10.9680 | 23.1530 | 7.0895 | 17.2990 | 15.3069 | 15.4760 | 15.4108 | 15.3179 | 15.2160
- | Combined | - | 11.3074 | - | - | 10.4607 | 10.4619 | 10.4606 | 10.4612 | 10.4621

We compare the generalization abilities of model trained on raw data, site-
removed data generated by ComBat (johnson2007adjusting;
pomponio2020harmonization), and site-independent data generated by flow-based
SCM for age prediction task as shown in Table 4. All models are firstly
trained on BLSA-1.5T (reference site) and then tested on BLSA-3T, UKBB, and
SHIP separately. The baseline column show validation MAEs directly trained on
each target sites. We find that model trained on the reference site with raw
data couldn’t generalize on the target sites. Models trained with site-removed
data generated by ComBat generalize much better compared to the one trained on
raw data. The performances of the variants (linear and GAM) of ComBat are
similar, except GAM shows lower MAE than linear when testing on BLSA-3T. The
model trained on site-independent data generated by Affine flow model has
slightly larger testing MAEs compared to ComBat models, though it couldn’t
capture the feature density as shown in Figure 3. Linear and quadratic Spline
flow models show relatively better performance than ComBat methods in all
target sites. For example, the linear Spline model achieves $5.9314$ MAE on
BLSA-3T, which is lower than $6.6203$ of ComBat-GAM. The quadratic Spline
model achieves $6.4992$ MAE and $15.2160$ MAE on UKBB and SHIP separately,
which is lower than $6.5294$ and $15.3069$ of ComBat-Linear respectively.
Finally, both ComBat and flow-based SCM show better prediction performance
when training on a combination of all studies in iSTAGING compared to the
original data.

###  5.6 Alzheimer’s Disease (AD) classification

Table 5: AD classification accuracy comparison in ADNI dataset. Baseline show results trained on target site directly. The raw, ComBat, and counterfactuals columns use the raw features, features generated by the ComBat algorithms, and features generated by the proposed method. | Study | TarOnly | SrcOnly | Transfer | IRM | ComBat | Flow-based SCM (ours)
---|---|---|---|---|---|---|---
Linear | GAM | Affine | L-Spline | Q-Spline
Source | ADNI-1 | - | 0.7355 | - | 0.7355 | 0.7273 | 0.7273 | 0.7521 | 0.7438 | 0.7355
Target | ADNI-2 | 0.7677 | 0.6968 | 0.7935 | 0.6452 | 0.7273 | 0.6942 | 0.6694 | 0.7521 | 0.7521
Source | ADNI-2 | - | 0.7677 | - | 0.7290 | 0.7851 | 0.7769 | 0.7603 | 0.7769 | 0.7934
Target | ADNI-1 | 0.7355 | 0.7190 | 0.7603 | 0.6777 | 0.6694 | 0.6529 | 0.6694 | 0.7355 | 0.7190
- | Combined | - | 0.8182 | - | - | 0.8141 | 0.7355 | 0.8099 | 0.8306 | 0.8430

Similar to age prediction, we evaluate model generalizibility when training on
the original data, site-removed data, and site-independent data from ADNI
dataset for AD classification task. Again, all models are firstly trained on
the source sites and then tested on the target sites respectively. The
baseline column show validation classification accuracies directly trained on
each target sites. We also find that there are significant generalization gaps
between source and target sites when training with raw data. Model trained
with site-removed data generated by ComBat-Linear show better testing accuracy
on ADNI-2 compared to the one trained on raw data, whereas it doesn’t perform
well in the other way (train on ADNI-2 and test on ADNI-1). ComBat-GAM models
show worse generalization ability compared to the ones train on raw data.
Similarly, models trained on site-independent data generated by Affine flow-
based SCM fail in improving the testing accuricies. We believe that the
feature location-scale correction process in ComBat methods removes the
disease information in raw data, which leads to unsatisfying performance. On
the other hand, the Affine flow model couldn’t capture the feature
distribution of the original data, which results in undesired generalization
ability. In contrast, both linear and quadratic Spline flow models show
significant testing accuracy improvement, which are close to the ones directly
trained on target sites. For example, both the linear and quadratic Spline
models achieve $0.7521$ testing accuracy on ADNI-2, which is close to $0.7677$
of baseline. The linear Spline model achieves $0.7355$ testing accuracy on
ADNI-1, which is exactly the same as the baseline model trained on ADNI-1.
Finally, the Spline flow-based SCM models show better generalization ability
when training on a combination of two phases of ADNI compared to the raw data.

###  5.7 Additional experiments using benchmark datasets

Table 6: Digit classification accuracy comparison in synthetic dataset. | Group | TarOnly | SrcOnly | Transfer | IRM | Flow-based SCM (ours)
---|---|---|---|---|---|---
Source | Plain | - | 0.9890 | - | 0.9820 | 0.9870
Target | Thin | 0.9790 | 0.9320 | 0.9970 | 0.9200 | 0.9230
Target | Thick | 0.9750 | 0.9160 | 0.9940 | 0.8980 | 0.9190
Source | Thin | - | 0.9790 | - | 0.9790 | 0.9840
Target | Plain | 0.9890 | 0.9220 | 0.9980 | 0.9470 | 0.9580
Target | Thick | 0.9750 | 0.5150 | 0.9940 | 0.6220 | 0.7020
Source | Thick | - | 0.9750 | - | 0.9120 | 0.9760
Target | Plain | 0.9890 | 0.9660 | 0.9980 | 0.8780 | 0.9640
Target | Thin | 0.9790 | 0.7530 | 0.9970 | 0.6380 | 0.7640

Digit classification with confounding variables.

thickness

horizontal/vertical

##  6 Discussion

In this paper, we tackle the confounding variations problem from a causal
view, showing that the vulnerability of discriminative learning algorithms can
be explained by the lack of causal reasoning. By explicitly modeling the
causal relationship of confoundings such as sex, age, and site, and imaging
data using our proposed flow-based SCM, we can generate counfouding-
independent data from counterfactual inference. Our model has demonstrated
improved robustness in both regression and classification tasks on a wide
range of real-world datasets. Furthermore, our proposed framework can be
enhanced by introducing more confoundings such as race and genetic
information. Finally, the framework can potentially benefit research in
explainable models by using counterfactual data.