# First-principles thermal equation of state of fcc iridium

Kai Luo<EMAIL_ADDRESS>Institute of Ultrafast Optical Physics, Department
of Applied Physics, Nanjing University of Science and Technology, Nanjing
210094, China  R. E. Cohen  Extreme Materials Initiative, Earth and Planets
Laboratory, Carnegie Institution for Science, 5241 Broad Branch Rd., N.W.,
Washington, DC 20015, USA  Ruifeng Lu  Institute of Ultrafast Optical Physics,
Department of Applied Physics, Nanjing University of Science and Technology,
Nanjing 210094, China

(Modified on Jan 5, 2023)

###### Abstract

The thermal equation of states for fcc iridium (Ir) is obtained from first-
principles molecular dynamics up to 3000 K and 540 GPa. The equation of state
(EoS) is globally fitted to a simplified free energy model and various
parameters are derived. The theoretical principal Hugoniot is compared with
shockwave experiments, where discrepancy suggests formation of new Ir phases.
A few representative EoS parameters, such as bulk modulus $K_{T}$, thermal
expansivity $\alpha$, Grüneisen parameter $\gamma$, and constant pressure
capacity $C_{P}$, Debye temperature, ${\Theta_{\mathrm{D}}}$ are computed to
compare with experimental data.

##  I Introduction

Iridium (Ir) is a 5d transition metal of the platinum group. It is the second-
densest metal with a density of 22.56 g/cm3 at ambient condition, only
slightly lower by about 0.12% than the densest metal osmium (Os). It has the
largest shear modulus, G=210 GPa among the face-centered cubic (fcc) metals.
The solid Ir remains in the fcc structure up to the melting point of 2719
K[1]. Due to its prominent thermophysical and mechanical properties and high
corrosion resistance, it is used in many technological applications, such as
crucibles, thermocouples, spark plugs, aircraft engine parts, and deep water
pipes. The lack of phase transitions, simple fcc structure, high melting
temperature, and non-reactivity, makes it ideal for experiments as a heater,
absorber, or standard for example in diamond-anvil cell (DAC) experiments, and
ideal for studying effects of compression on noble metals. Our understanding
of the properties of Ir is still limited, and fundamental research on it
remains of great interest.

With the advances in lab technologies, extreme conditions ($P>200$ GPa,
$T>2000$ K) become more and more amenable to study. Fundamental to all studies
at extreme conditions is the equation of state (EoS) that relates $P,V,T$, and
$U$ or $F$, where symbols of $P,V,T,U$, and $F$ stand for pressure, volume,
temperature, internal energy, and Helmholtz free energy. The earliest
investigation of iridium EoS dates back to 1937 by P.W. Bridgman up to 7
GPa[2, 3], followed by work of Schock and Johnson [4] and then of Akella[5] up
to 30 GPa. Cerenius and Dubrovinsky [6] measured the compressibility of Ir
using DAC up to 65 GPa. Later, Cynn et al. found that Ir has the second-lowest
compressibility of any element after Os from their DAC experiment up to 65
GPa, which was corroborated by first-principles calculations[7].

For the EoS diagrams, zero-temperature first-principles EoS can be
supplemented with finite-temperature vibrational entropies from the phonon
dispersions. Phonon frequencies can be calculated from finite differences, or
with the density-functional perturbation theory (DFPT) [8]. Thanks to the
development in the density functional theory toolkit, theoretical EoS for Ir
appeared in several experimental work [7, 9, 10, 11, 12]. However, these
theoretical EoS’s were limited to low temperatures (around 300 K) using static
calculations fitted to Birch-Murnaghan (BM) EoS [13]. Anharmonic lattice
vibrations were considered in Ref.[9], but the focus was the phase diagram and
phase stability. Anzellini et al. studied Ir up to 80 GPa and 3100K combining
in situ synchrotron X-ray diffraction using laser-heating DACs and density
functional theory calculations [14]. A comprehensive study covering a larger
range of temperatures and pressures has not been performed. Indeed, studying
other phases would be interesting, but in applications as a standard in
experiments, we focus on the fcc phase. In this work, we aim to provide the
EoS for fcc Ir up to 3000 K and 540 GPa in first-principles molecular dynamics
(FPMD).

##  II Theoretical EoS from FPMD

First principles methods have been widely adopted in the simulation of
condensed phases where no phenomenological parameters are needed. They give
access to a space of thermodynamic conditions, which are hard to reach for
experimental efforts and can be used to help calibrate experiments, where, for
example temperature data are not sometimes available at the desired
conditions. FPMD takes into account of anharmonic vibrations of ions directly
at finite temperatures through thermostatting. The electronic free energy is
given by the Mermin-Kohn-Sham density functional theory (DFT) [15, 16]. FPMD
becomes the most used tool for predicting the thermal EoS, subject to the
exchange-correlation free energy functional approximations[17, 18, 19].
Classical molecular dynamics is suitable for high temperatures above the Debye
temperature as it includes anharmonicity exactly, unlike other approaches. At
lower temperature deviations from the $P-V-T$ equation of state are small, but
heat capacities and high order properties such as thermal expansivity which
show strong quantum effects at temperature below the Debye temperature are
indeed less accurate.

###  II.1 FPMD details

First, we computed the static EoS of Ir at zero temperature. We used Quantum
Espresso ver. 6.7 throughout this work [20]. We used the scalar-relativisitic
(Garrity-Bennett-Rabe-Vanderbilt) GBRV ultrasoft pseudopotential [21] with
Perdew-Burke-Erzernhorf (PBE) exchange-correlation (xc) functional [22, 23].
The electronic configuration for the pseudopotential is
$[Xe]5p^{6.0}5d^{8.5}$. EoS was derived by fitting energy-volume curve in the
3rd-order BM equation. To validate the range of applicability of the
pseudopotential, we performed similar calculations in linearized augmented
planewave (LAPW) code Elk [24] and above two $P-V$ curves agree well up to 550
GPa (see Fig.1).

Figure 1:  Static equations of state from GBRV pseudopotential planewaves (QE)
and LAPW (elk) calculation are compared. The Vinet EoS was used to fit the
energy-volume curve. Inset figure shows the pressure difference and the
maximum is less than 10 GPa.

For the FPMD calculation, we prepared a cubic box containing 108 atoms in the
fcc structure. The energy cutoffs for planewaves and density are 80 Ry and 320
Ry, respectively. Energy is converged within 5 meV per atom. Only $\Gamma$
point was sampled. The bands are occupied according to the Fermi-Dirac
distribution at each temperature, and the number of bands are large enough to
guarantee the occupation number is smaller than $10^{-7}$ for the highest
occupied state. Early studies showed that the spin-orbit coupling does not
affect the EoS and hence we used spin-unpolarized DFT neglecting spin-orbit
coupling. Then conditions at a combination of lattice constants
$a/a_{0}=0.86,0.88,0.90,0.92,0.96,1.00$ ($a_{0}=3.801$ Å) and temperatures
$T=300,1000,1500,2000,2500,3000$ K were used in the simulations (see
conditions in Table 1). The time step is 20 a.u. (0.9676 fs). The equilibrated
time steps are more than 2000 to get the statistical means and standard
deviations, which give less than $1\%$ standard deviation. The ionic
temperature is regulated by the stochastic-velocity rescaling thermostat [25]
and no quantum corrections to the ionic motion are included.

Table 1: Pressure $P$ and internal energy per atom $U$ and its standard deviation of the means $P_{\mathrm{err}}$ and $U_{\mathrm{err}}$ are extracted from FPMD simulations of fcc Ir for a given temperature $T$ and atomic volume $V$ (or the mass density $\rho$). The global minimum of $U$ is underlined.  $T$ (K) |  $V$ (bohr3) |  $\rho$ (g/cm3) |  $P$ (GPa) |  $P_{\mathrm{err}}$ (GPa) |  $U$ (Ry) |  $U_{\mathrm{err}}$ (Ry)
---|---|---|---|---|---|---
300 | 61.01 | 35.31 | 527.6 | 0.01 | -181.353 | 0.00004
300 | 63.14 | 34.12 | 451.6 | 0.01 | -181.423 | 0.00004
300 | 67.54 | 31.89 | 326.1 | 0.02 | -181.538 | 0.00005
300 | 72.14 | 29.86 | 229.5 | 0.01 | -181.624 | 0.00004
300 | 81.97 | 26.28 | 99.3 | 0.02 | -181.729 | 0.00005
300 | 92.65 | 23.25 | 25.0 | 0.03 | -181.770 | 0.00006
1000 | 61.01 | 35.31 | 531.3 | 0.04 | -181.339 | 0.00012
1000 | 63.14 | 34.12 | 455.3 | 0.06 | -181.410 | 0.00021
1000 | 67.54 | 31.89 | 330.0 | 0.04 | -181.524 | 0.00014
1000 | 72.14 | 29.86 | 233.4 | 0.06 | -181.611 | 0.00017
1000 | 81.97 | 26.28 | 103.5 | 0.06 | -181.715 | 0.00016
1000 | 92.65 | 23.25 | 29.5 | 0.07 | -181.756 | 0.00018
1500 | 61.01 | 35.31 | 533.9 | 0.08 | -181.329 | 0.00039
1500 | 63.14 | 34.12 | 458.0 | 0.06 | -181.400 | 0.00022
1500 | 67.54 | 31.89 | 332.8 | 0.07 | -181.514 | 0.00021
1500 | 72.14 | 29.86 | 236.6 | 0.09 | -181.600 | 0.00033
1500 | 81.97 | 26.28 | 106.4 | 0.09 | -181.705 | 0.00027
1500 | 92.65 | 23.25 | 32.6 | 0.13 | -181.746 | 0.00034
2000 | 61.01 | 35.31 | 536.9 | 0.10 | -181.318 | 0.00037
2000 | 63.14 | 34.12 | 461.0 | 0.08 | -181.389 | 0.00027
2000 | 67.54 | 31.89 | 335.6 | 0.08 | -181.504 | 0.00027
2000 | 72.14 | 29.86 | 239.2 | 0.11 | -181.591 | 0.00034
2000 | 81.97 | 26.28 | 109.4 | 0.17 | -181.694 | 0.00054
2000 | 92.65 | 23.25 | 35.6 | 0.22 | -181.735 | 0.00076
2500 | 61.01 | 35.31 | 539.5 | 0.14 | -181.308 | 0.00038
2500 | 63.14 | 34.12 | 463.7 | 0.14 | -181.379 | 0.00051
2500 | 67.54 | 31.89 | 338.6 | 0.10 | -181.493 | 0.00032
2500 | 72.14 | 29.86 | 242.3 | 0.17 | -181.579 | 0.00058
2500 | 81.97 | 26.28 | 112.5 | 0.15 | -181.683 | 0.00046
2500 | 92.65 | 23.25 | 38.7 | 0.19 | -181.725 | 0.00053
3000 | 61.01 | 35.31 | 542.5 | 0.11 | -181.298 | 0.00037
3000 | 63.14 | 34.12 | 466.6 | 0.20 | -181.368 | 0.00072
3000 | 67.54 | 31.89 | 341.1 | 0.16 | -181.484 | 0.00057
3000 | 72.14 | 29.86 | 245.6 | 0.22 | -181.568 | 0.00075
3000 | 81.97 | 26.28 | 115.1 | 0.23 | -181.674 | 0.00068
3000 | 92.65 | 23.25 | 41.7 | 0.45 | -181.714 | 0.00141

###  II.2 Free energy model

We fit the Helmholtz free energy ($F$) as a function of $V$ and $T$, $F(V,T)$.
In FPMD, we have direct access to the variables of volume ($V$), temperature
($T$), pressure ($P$), and internal energy ($U$). Cohen and Gülseren [26]
studied the thermal EoS of tantalum (Ta) in full potential LAPW and mixed-
basis pseudopotential methods. An accurate high-temperature global EoS was
formed from the $T=0$ K Vinet isotherm and the thermal free-energy was fitted
by the polynomial expansion in $V$ and $T$ (see Eq. (11) in Ref. [26]). de
Koker and Stixrude [27] computed the free energy of MgO periclase and MgSiO3
perovskite using FPMD, where the excess free energy was fitted in a similar
expansion. Incorporating the Debye model [28], the total free energy is
approximated by the polynomial expansion up to order $N_{i},N_{j}$,

$F(V,T)=\sum_{\begin{subarray}{c}i,j=0\end{subarray}}^{N_{i},N_{j}}A_{ij}T^{i}(V^{-\frac{2}{3}})^{j}+F_{0}\,.$
(1)

Neglecting the zero-point motion,
$F_{0}=k_{B}T\left[-D_{3}(x)+3\ln(1-e^{-x})\right]$ where a dimensionless
parameter $x=\frac{{\Theta_{\mathrm{D}}}}{T}$ with Debye temperature
${\Theta_{\mathrm{D}}}$. $k_{B}$ is the Boltzmann constant. $D_{3}(x)$ is the
third order Debye function (see Appendix B). $A_{ij}$ are fitting coefficients
yet to be determined. For comparison, we mention the classical model, where
$F_{0}=-3k_{B}T\ln T$, with $T\ln T$ giving the proper classical behavior at
low temperatures. That is $C_{V}=3k_{B}$ and $S=-\infty$ at 0 K. The Debye
temperature ${\Theta_{\mathrm{D}}}$ cannot be determined from the
$U(T,V),P(T,V)$ data from the classical molecular dynamics, so we obtain
${\Theta_{\mathrm{D}}}$ from the RMS displacements (see below). For
simplicity, the Debye temperature at $P=0$ GPa, $T$=300 K, is used.

##  III Results

We obtained the equilibrated quantities from FPMD, where $U,P$ includes the
ionic kinetic energy and ideal gas pressure, respectively. We subtracted each
internal energy by the global minimum, as only the energy difference matters.
The pressure and internal energy are $P=-\left(\frac{\partial F}{\partial
V}\right)_{T},U=F+TS=F-T\left(\frac{\partial F}{\partial T}\right)_{V}$.
$(U,P)$ data are grouped as a pair and fitted together to avoid bias between
these two quantities. The fitting was performed using the weighted least-
square fit with the lm function including offset in R language. Internal
energy $U$ and pressure $P$ were fitted simultaneously to $F(V,T)$.
$w=1/\Delta^{2}$ is set for the weight, where $\Delta$ is the standard
deviation of $U$ and $P$. We fitted Eq. (1) with $N_{i}=2,N_{j}=3$. We
analyzed the MD trajectories using the code VMD, and computed the root-mean-
square displacement (RMSD) for each run. From this we can obtain the effective
Debye temperature ${\Theta_{\mathrm{D}}}$ using:

$\langle
u^{2}\rangle=\frac{3h^{2}}{4\pi^{2}Mk_{B}{\Theta_{\mathrm{D}}}}\left(\frac{D_{1}({\Theta_{\mathrm{D}}}/T)}{{\Theta_{\mathrm{D}}}/T}+\frac{1}{4}\right)\,,$
(2)

where $\vec{u},M,h$ are the displacement vector, the ion mass, the Planck
constant, and $D_{1}$ is the first order Debye function. The quantum
correction term $\frac{1}{4}$ shall be omitted in the classical treatment.
Since the phonon density of states is not exactly Debye-like, this is the
effective Debye temperature for the second moment of the vibrational density
of states (VDOS) , not exactly equal to the thermodynamic Debye temperature
[29].

Figure 2:  The residuals of the fit to Eq. (1) for the EoS of fcc iridium are
shown in a) and b). $\Delta U=U-U_{min}$, where $U_{min}$ is the minimum in
the dataset underlined in Table 1. The fitted curves are compared against the
dataset in c) and d). $\Delta U$ and thermal pressure $P-P_{\mathrm{300K}}$
from the fit align well against the dataset.

The residual is the deviation between the target function and the sample mean.
From Fig. 2, we observe that the residuals are randomly distributed across the
volume range. The absolute value of residuals for $U$ and $P$ are less than
0.002 Ry and 1.0 GPa (except for data point at 3000 K). For the internal
energy, a global minimum $U_{min}$ is subtracted from the dataset. On the
scale of half Ry, the internal energy is well represented. As for the
pressure, we computed the pressure differences with respect to the $T=300$ K
reference and the fitted curves aligned with the dataset. Only the $T=3000$ K
fit is slightly off. The resultant fitting coefficients in atomic unit for
both the Debye model and the classical model are tabulated (see Table 2). The
statistical summary from lm function is included in the Appendix A (see Fig.
12).

Table 2: Coefficient matrix $A$ in atomic units for the choice of $F_{0}$, the Debye model and the classical model. The root mean squared errors (RMS) in the fitting for the pressure $P$ (in GPa) and the internal energy $U$ (in mRy) are listed.  | $A_{ij}$ |  $P$ RMS (GPa) |  $U$ RMS (mRy)
---|---|---|---
Debye model |  1.904 -62.84 88.81 8175 0 0.008698 -0.1175 0.5574 $6.681\times 10^{-09}$ $-3.752\times 10^{-07}$ $6.321\times 10^{-06}$ $-3.439\times 10^{-05}$ | 0.5380 | 0.706
Classical model |  1.902 -62.83 88.68 8176 0 0.008676 -0.1171 0.5549 $6.567\times 10^{-09}$ $-3.667\times 10^{-07}$ $6.166\times 10^{-06}$ $-3.347\times 10^{-05}$ | 0.5378 | 0.661

###  III.1 $P-V-T$ EoS

The equilibrium atomic volume ($P=0$ GPa) at 300 K is 14.559 Å3, 2.9% larger
than the experimental value 14.145 Å3. The overestimation of the lattice
constants is expected for the PBE exchange-correlation functional.
Experimental $P-V$ curves of 300 K isotherm are readily compared with our
theoretical predictions. Pressures measured by Akella et al. [5] are
underestimated for compression (see Fig. 3), $\Delta V/V_{0}$ larger than 0.05
with $\Delta V=V_{0}-V$. Overall the theoretical 300 K isotherm agrees well
with the experiments within the uncertainty especially when the compression is
smaller than 0.15 ($P<70$ GPa) [6, 30, 10]. In contrast, the 3rd order BM fit
done by Monteseguro et al. [10] sits along our 1000 K isotherm for compression
$>0.15$, and reflects the inadequecy of BM EoS at high compression. For
comparison, we have also included the FPMD and experimental study of Anzellini
et al. [14]. Their $P-V-T$ curves (both theory and experiments) below 80 GPa
are obtained using the EoSFit7 package with ingredients such as the third-
order BM EoS for the isothermal part. Their FPMD used the local density
approximations and smaller energy cutoff (300 eV). Isotherm of 0 K compared
well against our 300 K curve at low compression but not at high compression
(compression $>$ 0.90). Similar for the isotherms of 1000 K and 3000 K. The
shock-wave experiment by Al’tshuler et al. [31, 32, 33] exhibits quite
distinct behavior in the $P-V$ curve. Around compression of 0.1, the
temperature is pinned slightly above the isotherm of 1000 K and at compression
of 0.22 the temperature is close to the 3000 K isotherm. The high compression
pressure ($\approx 600$ GPa) of Al’tshuler et al. was mistakenly reported in
Ref. [10]. The recent shockwave experimental work by Khishchenko [34] is also
compared. We observe the room temperature isotherm of recent work by
Khishchenko et al. align almost perfectly with our EoS data. The data by
Monteseguro et al. runs along our 1000 K isotherm for compression over 0.1. It
is well-known that dynamic compression experiment lead to a temperature rise.
Contrary to the claim that the temperature effect is negligible by Monteseguro
et al., [10] we believe the temperatures increase (not measured) along the
shock compression $P-V$ curve is significant from our predicted EoS.

Figure 3:  Theoretical and experimental EoS’s of fcc iridium are compared. The
DAC experimental data for Yusenko [30] Monteseguro[10], Cerenius[6], and
Akella[5] were compared at 300 K, where the BM EoS’s were available.
Experimental and theoretical data from Anzellini [14] are included (pink,
green, and yellow lines). The shockwave data (red cross) of Al’tshuler were
taken from Ref. [10]. Room temperature isotherm of Khishchenko shock
experiments [34] agrees very well against our FPMD results.  Figure 4:  The
thermal pressure of fcc iridium is roughly linear in $T$.

Thermal pressure measures the pressure change upon temperature increase at
constant volume, $P_{th}(V,T)=P(V,T)-P(V,T_{0})$. The thermal pressure is
quite linear in $T$ given that $\alpha K_{T}$ ($\alpha$ and $K_{T}$ are the
thermal expansivity and the bulk modulus) is constant in the classical regime
(above the Debye temperature), expressed as

$P_{th}(V,T)=\int_{T_{0}}^{T}\\!\\!dT\left(\frac{\partial P}{\partial
T}\right)_{V}=\int_{T_{0}}^{T}\\!\\!dT\,\alpha K_{T}\,.$ (3)

An oversimplified linear equation (see Fig.4) could be given to the thermal
pressure $P_{th}(T)=\lambda T$, with $\lambda=0.0056$ GPa/K for the
equilibrium volume. One could also see the volume dependence is weak from the
bottom right panel of Fig. 2.

The equilibrium bulk modulus $B_{0}$ (or inverse compressibility at room
temperature and zero pressure) is an important parameter in the EoS formula,
such as the Vinet EoS [35, 36]. The fitted bulk modulus is compared against
earlier studies (see Table. 3). We note that Cerenius and Dubrovinsky [6]
obtained similar bulk modulus, $354$ GPa versus $306$ GPa, by fitting the
second order BM EoS with constraint $B_{0}^{\prime}=4$, or third-order BM EoS
without constraint both using experimental equilibrium volume. Park et al.
obtained the bulk modulus of 399 GPa and 344 GPa for the LDA and GGA
functional in DFT, respectively [37]. We note that $B_{0}$ from our fit is
close to the accepted value of about 365 GPa and evidently smaller than Cynn’s
value 383 GPa [7]. The parameter $B_{0}^{\prime}$ from our model is $5.3$.

Table 3: Experimental equilibrium volume $V_{0}$ (Å3 per atom), bulk modulus $B_{0}$ (GPa), and $B_{0}^{\prime}$ at room temperature are compared against reported theoretical results. Data and method are briefly summarized, and the original references are given.  Method description | $V_{0}$ | $B_{0}$ | $B_{0}^{\prime}$ | References
---|---|---|---|---
| (Å3/at) | (GPa) |  |
Exp. data fitted to 3rd-order BM EoS | 14.120 | 339 | 5.3 | Monteseguro et al. [10]
Exp. data fitted to 3rd-order BM EoS | 14.145 | 383 | 3.1 | Cynn et al. [7]
Exp. data fitted to |  |  |  | Cerenius and Dubrovinsky[6]
2nd-order BM EoS, with $B_{0}^{\prime}=4$ | 14.173 (exp. value) | 354 | 4.0 |
3rd-order BM EoS, without constraint | 14.173 (exp. value) | 306 | 6.8 |
DFT data fitted to BM EoS |  |  |  | Park et al. [37], Table 1 and 2
PAW LDA | 13.925 | 399 |  |
PAW GGA | 14.524 | 344 |  |
FPMD data fitted to 3rd-order BM EoS | 14.150 | 366 | 5.0 | Burakovsky et al. [9]
FPMD data fitted to our EoS | 14.559 | 361 | 5.3 | This work

###  III.2 Shock compression

High pressure high temperature conditions are generated by laser heating [38]
or resistive heating [39] in a DAC or by laser or gas gun [40] driven shock
compression. Strong shocks obey the Rankine-Hugoniot,

$U-U_{0}+\frac{1}{2}(P+P_{0})(V-V_{0})=0\,.$ (4)

Since the analytical expression for $U,P$ as a function of $V,T$ is known, for
each volume $V$, we solve Eq. (4) by searching its root $T$ given the
experimental value $V_{0},T_{0}$.

Figure 5:  The principle Hugoniot curve from our theoretical EoS (solid black)
of fcc Ir is compared against the shockwave experimental data (symbols). Round
green is for the STAR Hugoniot [41], square blue for the LANL Hugoniot [42],
and diamond orange for Al’tshuler [31, 32], respectively. Computed temperature
from the Eq. (4) is shown in red. The dashed extension is beyond the
simulation domain.

We compared our predicted principle Hugoniot with available shock experimental
data from several facilities [43, 41] (Fig. 5). Our theoretical principle
Hugoniot agrees well with that from earlier data of Al’tshuler and LANL March,
as well as more recent data of STAR Hugoniot, for $P<200$ GPa. Above 200 GPa,
our predicted pressure is higher than that of LANL and STAR but lower than
Al’tshuler’s. Our theoretical Hugoniot below 3000 K (shock temperature) are
fairly reliable which correspond to pressure less than 200 GPa. The shock
temperature, calculated as the solution to Eq. (4), is shown.

###  III.3 Equation of state parameters

Thermal EoS parameters such as thermal expansivity $\alpha$, isothermal
compressibility $\beta_{T}$, Grüneisen parameter $\gamma$, and the heat
capacity $C_{V}$ and $C_{P}$ are obtained by differentiation and algebraic
manipulation of Eq. (1). We now discuss some of these parameters.

Figure 6:  Parameter $\alpha K_{T}$ as a function of temperature for various
pressures.

As expected from the thermal pressure, $\alpha K_{T}$ is weakly dependent on
the volume and temperature (see Fig. 6). The Grüneisen parameter

$\gamma=V\left(\frac{\partial P}{\partial U}\right)_{V}=V\frac{\alpha
K_{T}}{C_{V}}$ (5)

is another important parameter. It is used in the Mie-Grüneisen EoS, where
$\gamma$ is assumed independent of temperature. The span of $\gamma$ as a
function of temperature reduces when the pressure increases (see Fig. 7). For
high compressions, it is indeed fairly temperature independent.

Figure 7:  Grüneisen parameter $\gamma$ as a function of temperature for
various pressures.  Figure 8:  The volumetric thermal expansivity $\alpha$ of
fcc iridium as a function of pressure at various temperatures. Halvorson and
Wimber measured the linear thermal expansion as
$\alpha_{t}=a_{0}+a_{1}t+a_{2}t^{2}+a_{3}t^{3}$ with $a_{0}=6.167\times
10^{-6},a_{1}=3.038\times 10^{-9},a_{2}=-0.8448\times
10^{-12},a_{3}=0.5852\times 10^{-15}$, for $t$ expressed in ∘C [44] at ambient
pressure (see inset solid line), where one can show
$\alpha=3(L_{0}/L_{t})\alpha_{t}$ for isotropic materials with reference
length $L_{0}$.

The volumetric thermal expansivity $\alpha=-1/V(\partial V/\partial T)_{P}$
for isotropic materials is three times the linear thermal expansivity
coefficient $\alpha_{L}$, $\alpha=3\alpha_{L}$. $\alpha$ in Fig. 8 is
essentially temperature independent but rather volume sensitive. Our
theoretical prediction is below the reported experimental value [44], but it
is noted that around room temperature, our theory prediction gives the right
thermal expansion coefficient. The expansivity has downsized by a factor of 4
when the pressure goes to 300 GPa.

Figure 9:  Constant pressure heat capacity $C_{P}$ as a function of
temperature (top panel) and pressure (bottom panel). Experimental value at
ambient condition is 25.10
$\mathrm{J}\text{\,}{\mathrm{K}}^{-1}\text{\,}{\mathrm{mol}}^{-1}$[45]. The
classical model (blue dash-dotted) for 0 GPa starts to deviate for $T$ below
500 K, and approaches the classical limit as $T\to 0$ K. Inset shows the
constant volume heat capacity $C_{V}$. Heat capacity largely reduces when
pressure goes up.

The heat capacity at high pressures are almost linear above 500 K, see Fig. 9.
The higher the temperature, the slope of $C_{P}$ is smaller. At 0 GPa, the
predicted value 25.68
$\mathrm{J}\text{\,}{\mathrm{K}}^{-1}\text{\,}{\mathrm{mol}}^{-1}$is fairly
accurate and only 2.3% larger than the experimental heat capacity 25.10
$\mathrm{J}\text{\,}{\mathrm{K}}^{-1}\text{\,}{\mathrm{mol}}^{-1}$[45], given
that we used the formula, $C_{P}=C_{V}(1+T\alpha\gamma)$ (see Appendix A) with
errors in $C_{V},\alpha,$ and $\gamma$.

The RMSD is a critical quantity for the analysis of the phonon vibrations.
Moseley et al. [46] presented temperature-dependent inelastic neutron
scattering (INS) experiments as well as quasi-harmonic density functional
theory calculations to study the thermodynamic properties of Ir. Our FPMD RMSD
at 300 K agrees particularly well with their experimental findings. Their
reported $\langle u^{2}\rangle$ at higher temperatures (T=673 K and 823 K, see
Table I of Ref. [46]) however, are higher than our FPMD predictions. This is
reasonable since our NVT ensembles at these temperatures lead to higher
pressure and confined vibrations. It is worth to note that their phonon
density of states (PDOS) integrates to 1, and is not fitted well at higher
energies. Further investigation of PDOS with a FPMD simulation to compare
against the experiments may give insight to the anharmonic effects. Debye
temperatures of isochores using Eq. (2) exhibit weak temperature
dependence(see Fig. 11).

Figure 10:  RMSD as a function of temperature. Experimental data was obtained
by Moseley et al. [46] for constant pressure.  Figure 11:  Debye temperature
${\Theta_{\mathrm{D}}}$ as a function of temperature. Zero point motion is not
included to obtain the Debye temperature. In the high temperature, classical
region, our RMSD are classical from classical FPMD (Fig. 10) but when an
effective classical Debye temperature is derived to model the RMSD, there is a
large change with decreasing temperature into to quantum regime in the Debye
model. A 4th order polynomial fit is done only for higher temperature due to
the classical treatment to the ions.

##  IV Conclusions

We have performed a series of FPMD simulations for the fcc Ir at conditions up
to 3000 K and 540 GPa. By using a simplified model for the free-energy as a
function of temperature and volume and the statistical average quantities
internal energy and pressure $(U,P)$, the thermal EoS is obtained by globally
fitting to the model. We have compared previous experimental EoS’s and
provided the thermal EoS up to 3000 K, and 540 GPa. The $P-V-T$ curve is
reasonably agreeing with the fitted BM EoS at low compression but differs at
high compression. Our first-principles EoS accords with the most recent
shockwave experiment work by Khishchenko We find that $\alpha K_{T}$ and the
thermal pressure is quite constant from its dependence in temperature and
which turns out to be true for a wide range of materials. We have shown some
representative derived thermal parameters against available experiments and
found agreements and discrepancies. Further work might resolve these
discrepancies.

##  V Acknowledgments

The work was done under the auspices of the US National Science Foundation
CSEDI grant EAR-1901813 to R.E.C. and the National Natural Science Foundation
of China (Grant No. 12104230) to K.L.; R.E.C. is supported by the Carnegie
Institution for Science and gratefully acknowledges the Gauss Centre for
Supercomputing e.V. for funding this project by providing computing time on
the GCS Supercomputer Supermuc-NG at Leibniz Supercomputing Centre. All the
FPMD calculations were performed on Supermuc-NG.

##  Appendix A Thermodynamic relations

Once the Helmholtz free energy $F=F(V,T)$ is known as a function of volume
($V$) and temperature ($T$), the following thermodynamical quantities can be
obtained from it [47]:

$\displaystyle P$ $\displaystyle=$ $\displaystyle-\left(\frac{\partial
F}{\partial V}\right)_{T}\,,$ (6) $\displaystyle S$ $\displaystyle=$
$\displaystyle-\left(\frac{\partial F}{\partial T}\right)_{V}\,,$ (7)
$\displaystyle K_{T}$ $\displaystyle=$
$\displaystyle\beta_{T}^{-1}=-V\left(\frac{\partial P}{\partial
V}\right)_{T}=V\left(\frac{\partial^{2}F}{\partial V^{2}}\right)_{T}\,,$ (8)
$\displaystyle C_{V}$ $\displaystyle=$ $\displaystyle T\left(\frac{\partial
S}{\partial T}\right)_{V}=-T\left(\frac{\partial^{2}F}{\partial
T^{2}}\right)_{V}\,,$ (9) $\displaystyle\alpha K_{T}$ $\displaystyle=$
$\displaystyle-\left(\frac{\partial^{2}F}{\partial T\partial V}\right)\,,$
(10) $\displaystyle\gamma$ $\displaystyle=$ $\displaystyle
V\left(\frac{\partial P}{\partial U}\right)_{V}=V\frac{\alpha
K_{T}}{C_{V}}\,,$ (11) $\displaystyle\frac{C_{P}}{C_{V}}$ $\displaystyle=$
$\displaystyle\frac{K_{S}}{K_{T}}=1+T\alpha\gamma\,,$ (12) $\displaystyle U$
$\displaystyle=$ $\displaystyle F+TS\,,$ (13)

where pressure, entropy, isothermal compressibility (its inverse is the bulk
modulus $K_{T}$), constant volume molar heat capacity, volumetric expansion
coefficient, Grüneisen parameter, internal energy are denoted by
$P,S,\beta_{T},C_{V},\alpha,\gamma,U$. $C_{P}$ is the constant pressure molar
heat capacity. The parameter $B_{0}^{\prime}$ (or
$K_{0}^{\prime}=\left.\frac{\partial K}{\partial P}\right|_{P=0}$) can thus be
computed using above relations:

$\displaystyle K^{\prime}$ $\displaystyle=$ $\displaystyle\frac{\partial
K}{\partial P}$ (14) $\displaystyle=$ $\displaystyle\left(\frac{\partial
K}{\partial T}\right)_{V}\left(\frac{\partial T}{\partial
P}\right)_{V}+\left(\frac{\partial K}{\partial
V}\right)_{T}\left(\frac{\partial V}{\partial P}\right)_{T}$ $\displaystyle=$
$\displaystyle\left(\frac{\partial K}{\partial
T}\right)_{V}\frac{1}{\left(\frac{\partial P}{\partial
T}\right)_{V}}+\left(\frac{\partial K}{\partial
V}\right)_{T}\frac{1}{\left(\frac{\partial P}{\partial V}\right)_{T}}$
$\displaystyle=$ $\displaystyle\left(\frac{\partial K}{\partial
T}\right)_{V}\frac{1}{\left(\frac{\partial P}{\partial
T}\right)_{V}}-\left(\frac{\partial K}{\partial
V}\right)_{T}\frac{V}{K_{T}}\,.$

Figure 12:  The fit summary plot from the lm function in R.

##  Appendix B Debye Model

The Helmholtz free energy $F$ of a vibrating lattice at volume $V$ and
temperature $T$, can be approximated as

$F(V,T)=E(V)+F_{\mathrm{vib}}(V,T)+F_{el}(V,T)\,,$ (15)

where $F_{\mathrm{vib}}$ is the vibrating energy of the lattice and $F_{el}$
is the thermal electronic free energy which is typically negligible. Moruzzi
et al. [28] proposed an empirical Debye model with

$F_{\mathrm{vib}}=k_{B}T\left[-D_{3}(x)+3\ln(1-e^{-x})\right]+\frac{9}{8}k_{B}{\Theta_{\mathrm{D}}}\,,$
(16)

with Debye temperature ${\Theta_{\mathrm{D}}}$ and dimensionless parameter
$x=\frac{{\Theta_{\mathrm{D}}}}{T}$. The last term is the zero-point energy.
$D_{3}(x)$ is the third order Debye function. The $n$th order Debye function
is defined as

$D_{n}(x)=\int_{0}^{x}\frac{t^{n}}{e^{t}-1}dt,x\geq 0\,,$ (17)

where $n$, a non-negative integer, is the order of the Debye function. The
vibrational entropy is

$S_{\mathrm{vib}}=4k_{B}D(x)-3k_{B}\ln(1-e^{-x})\,.$ (18)

Neglecting the zero-point motion, the vibrational internal energy
$U_{\mathrm{vib}}$ thus can be obtained

$U_{\mathrm{vib}}=F_{\mathrm{vib}}+TS_{\mathrm{vib}}=3k_{B}TD(x)\,.$ (19)

## References

  * Arblaster [2010] J. W. Arblaster, Platinum Metals Review 54, 93 (2010).
  * Bridgman [1937] P. W. Bridgman, Proceedings of the American Academy of Arts and Sciences 71, 387 (1937).
  * Bridgman [1952] P. W. Bridgman, Proceedings of the American Academy of Arts and Sciences 81, 165 (1952).
  * Schock and Johnson [1971] B. Schock and Q. Johnson, Fiz . Metal. Metalloved. 31, 1101 (1971).
  * Akella [1982] J. Akella, Journal of Physics and Chemistry of Solids 43, 941 (1982).
  * Cerenius and Dubrovinsky [2000] Y. Cerenius and L. Dubrovinsky, Journal of Alloys and Compounds 306, 26 (2000).
  * Cynn  _et al._ [2002] H. Cynn, J. E. Klepeis, C.-S. Yoo, and D. A. Young, Phys. Rev. Lett. 88, 135701 (2002).
  * Baroni  _et al._ [2001] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  * Burakovsky  _et al._ [2016] L. Burakovsky, N. Burakovsky, M. J. Cawkwell, D. L. Preston, D. Errandonea, and S. I. Simak, Phys. Rev. B 94, 094112 (2016).
  * Monteseguro  _et al._ [2019] V. Monteseguro, J. Sans, V. Cuartero, F. Cova, I. A. Abrikosov, W. Olovsson, C. Popescu, S. Pascarelli, G. Garbarino, H. J. M. Jönsson, _et al._ , Scientific reports 9, 1 (2019).
  * Fang  _et al._ [2010] H. Fang, B. Liu, M. Gu, X. Liu, S. Huang, C. Ni, Z. Li, and R. Wang, Physica B: Condensed Matter 405, 732 (2010).
  * Kaptay [2015] G. Kaptay, J. Mater. Sci. 50, 678 (2015).
  * Birch [1947] F. Birch, Phys. Rev. 71, 809 (1947).
  * Anzellini  _et al._ [2021] S. Anzellini, L. Burakovsky, R. Turnbull, E. Bandiello, and D. Errandonea, Crystals 11 (2021), 10.3390/cryst11040452.
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Mermin [1965] N. D. Mermin, Phys. Rev. 137, A1441 (1965).
  * Karasiev  _et al._ [2014] V. V. Karasiev, T. Sjostrom, J. Dufty, and S. B. Trickey, Phys. Rev. Lett. 112, 076403 (2014).
  * Karasiev  _et al._ [2016] V. V. Karasiev, L. Calderín, and S. B. Trickey, Phys. Rev. E 93, 063207 (2016).
  * Karasiev  _et al._ [2018] V. V. Karasiev, J. W. Dufty, and S. B. Trickey, Phys. Rev. Lett. 120, 076401 (2018).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Garrity  _et al._ [2014] K. F. Garrity, J. W. Bennett, K. M. Rabe, and D. Vanderbilt, Computational Materials Science 81, 446 (2014).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
  * Perdew  _et al._ [1997] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 78, 1396 (1997).
  * [24] http://elk.sourceforge.net/, .
  * Bussi  _et al._ [2007] G. Bussi, D. Donadio, and M. Parrinello, The Journal of Chemical Physics 126, 014101 (2007).
  * Cohen and Gülseren [2001] R. E. Cohen and O. Gülseren, Phys. Rev. B 63, 224101 (2001).
  * De Koker and Stixrude [2009] N. De Koker and L. Stixrude, Geophysical Journal International 178, 162 (2009), https://academic.oup.com/gji/article-pdf/178/1/162/5992596/178-1-162.pdf .
  * Moruzzi  _et al._ [1988] V. L. Moruzzi, J. F. Janak, and K. Schwarz, Phys. Rev. B 37, 790 (1988).
  * Wallace [1965] D. C. Wallace, Phys. Rev. 139, A877 (1965).
  * Yusenko  _et al._ [2019] K. V. Yusenko, S. Khandarkhaeva, T. Fedotenko, A. Pakhomova, S. A. Gromilov, L. Dubrovinsky, and N. Dubrovinskaia, Journal of Alloys and Compounds 788, 212 (2019).
  * Al'tshuler and Bakanova [1969] L. V. Al'tshuler and A. A. Bakanova, Soviet Physics Uspekhi 11, 678 (1969).
  * Al’Tshuler  _et al._ [1981] L. Al’Tshuler, A. Bakanova, I. Dudoladov, E. Dynin, R. Trunin, and B. Chekin, Journal of Applied Mechanics and Technical Physics 22, 145 (1981).
  * Nemoshkalenko  _et al._ [1988] V. Nemoshkalenko, V. Y. Mil’Man, A. Zhalko-Titarenko, V. Antonov, and Y. L. Shitikov, Soviet Journal of Experimental and Theoretical Physics Letters 47, 295 (1988).
  * Khishchenko [2022] K. V. Khishchenko, Journal of Physics: Conference Series 2154, 012009 (2022).
  * Vinet  _et al._ [1986] P. Vinet, J. Ferrante, J. R. Smith, and J. H. Rose, Journal of Physics C: Solid State Physics 19, L467 (1986).
  * Vinet  _et al._ [1987] P. Vinet, J. R. Smith, J. Ferrante, and J. H. Rose, Phys. Rev. B 35, 1945 (1987).
  * Park  _et al._ [2015] J. Park, B. D. Yu, and S. Hong, Current Applied Physics 15, 885 (2015).
  * Meng  _et al._ [2006] Y. Meng, G. Shen, and H. K. Mao, J. Phys.: Condens. Matter 18, S1097 (2006).
  * Boehler [1993] R. Boehler, Nature 363, 534 (1993).
  * Mitchell and Nellis [1981] A. C. Mitchell and W. J. Nellis, Review of Scientific Instruments 52, 347 (1981), https://doi.org/10.1063/1.1136602 .
  * Seagle  _et al._ [2019] C. Seagle, W. Reinhart, S. Alexander, J. Brown, and J.-P. Davis, “Shock compression of iridium,” 21st Biennial Conference of the APS Topical Group on Shock Compression of Condensed Matter (2019).
  * Marsh [1980] S. P. Marsh,  _LASL shock Hugoniot data_ (University of California press, 1980).
  * Fortov  _et al._ [2013] V. Fortov, L. Altshuler, R. Trunin, and A. Funtikov,  _High-Pressure Shock Compression of Solids VII: Shock Waves and Extreme States of Matter_, Shock Wave and High Pressure Phenomena (Springer New York, 2013).
  * Halvorson and Wimber [1972] J. J. Halvorson and R. T. Wimber, Journal of Applied Physics 43, 2519 (1972), https://doi.org/10.1063/1.1661553 .
  * Lide [2004] D. R. Lide,  _CRC handbook of chemistry and physics_ , Vol. 85 (CRC Press, 2004).
  * Moseley  _et al._ [2020] D. H. Moseley, S. J. Thébaud, L. R. Lindsay, Y. Cheng, D. L. Abernathy, M. E. Manley, and R. P. Hermann, Phys. Rev. Materials 4, 113608 (2020).
  * Callen [1998] H. B. Callen,  _Thermodynamics and an Introduction to Thermostatistics_ (American Association of Physics Teachers, 1998).


# Regularized regression on compositional trees with application to MRI
analysis

Bingkai Wang1 Brian S. Caffo1 Xi Luo2 Chin-Fu Liu3 Andreia V. Faria4 Michael
I. Miller3 Yi Zhao5 and for the Alzheimer’s Disease Neuroimaging
Initiative111Data used in preparation of this article were obtained from the
Alzheimer’s Disease Neuroimaging Initiative (ADNI) database
(adni.loni.usc.edu). As such, the investigators within the ADNI contributed to
the design and implementation of ADNI and/or provided data but did not
participate in analysis or writing of this report. A complete list of ADNI
investigators can be found at: http://adni.loni.usc.edu/wp-
content/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf
1Department of Biostatistics, Johns Hopkins Bloomberg School of Public Health
2Department of Biostatistics and Data Science, The University of Texas Health
Science Center at Houston
3Center for Imaging Science, Biomedical Engineering, Johns Hopkins University
4Department of Radiology, Johns Hopkins University School of Medicine
5Department of Biostatistics, Indiana University School of Medicine

###### Abstract

A compositional tree refers to a tree structure on a set of random variables
where each random variable is a node and composition occurs at each non-leaf
node of the tree. As a generalization of compositional data, compositional
trees handle more complex relationships among random variables and appear in
many disciplines, such as brain imaging, genomics and finance. We consider the
problem of sparse regression on data that are associated with a compositional
tree and propose a transformation-free tree-based regularized regression
method for component selection. The regularization penalty is designed based
on the tree structure and encourages a sparse tree representation. We prove
that our proposed estimator for regression coefficients is both consistent and
model selection consistent. In the simulation study, our method shows higher
accuracy than competing methods under different scenarios. By analyzing a
brain imaging data set from studies of Alzheimer’s disease, our method
identifies meaningful associations between memory declination and volume of
brain regions that are consistent with current understanding.

Keywords: composition, hierarchical tree, regularized regression.

##  1 Introduction

Compositional data refer to a type of data where data points are non-negative
and the data vector of each subject or observational unit sums up to one.
Compositional data appear in many disciplines, such as econometrics (Mullahy,,
2015), geology (Pawlowsky-Glahn and Egozcue,, 2006) and epidemiology (Leite,,
2016).

In the area of brain imaging, structural magnetic resonance imaging (MRI) and
anatomical brain segmentation produce compositional data. For example, using
3-dimensional images acquired via structural MRI, a five-step brain
segmentation introduced by Mori et al., (2016) can partition the whole brain
into regions at five granularity levels. At the most coarse level, the whole
brain is segmented into telencephalon (left and right), diencephalon (left and
right), metencephalon, mesencephalon and cerebrospinal fluid (CSF). At the
finest level, the whole brain is segmented into 236 brain regions. The
compositional data are then the fractional volumes of the 236 brain regions
relative to the intracranial volume (ICV).

In addition to composition, the volumetric data have a tree structure. In the
first step of segmentation, the whole brain is partitioned into 7 brain
regions. In the second step, each of the 7 brain regions created by the first
step is further partitioned into smaller regions, which can be thought of as
tree branching. Applied to all brain segmentation steps, this analogy makes a
tree structure that is rooted at the whole brain and has 236 leaves, which are
the brain regions at the finest segmentation. The tree structure is shown in
Figure 1. A key feature of this tree structure is that the volume of a brain
region is equal to the combined volume of its subregions (after one
segmentation), which introduces extra composition among variables. We refer to
this data structure as “compositional tree”. We note that the structure of
compositional data is a special case of compositional trees, which only have
leaves and a root.

Figure 1: Compositional tree structure of the MRI data example. The tree is
rooted at the whole brain. A brain region with suffix “_L” or “_R” indicates
that the region is in the left or right hemisphere of the brain.

Compositional trees appear in many disciplines. For example, Wang and Zhao,
(2017) presented a compositional tree of microbiome data, where the
compositional tree is formed by bacterial taxa at multiple taxonomic levels.
Another example is the fractional market capitalization of stocks in the S&P
500 index (relative to the total market capitalization of S&P 500), where all
500 stocks are partitioned into 11 sectors and each sector is further broken
down into industries according to the Global Industry Classification Standard
(MSCI,, 2020). The fractional market capitalization of a sector (or industry)
is the summation of fractional market capitalization of stocks that are
categorized into this sector (or industry). Compared with compositional data,
compositional trees provide more information about the relationships among
variables and suggest grouping effects at different levels.

Although methods for analyzing compositional data or tree-structured data have
been developed, little is known about how to deal with compositional trees.
Lin et al., (2014); Fiksel et al., (2020); Ma and Zhang, (2020) studied
regression methods for compositional data with or without regularization, but
their results cannot be directly generalized to handle compositional trees.
Kim and Xing, (2012) proposed a tree lasso for estimating a sparse multi-
response regression function, which did not consider compositional data. To
the best of our knowledge, the primary competitive work is Wang and Zhao,
(2017), which developed a tree-guided regularization method for structured
sub-composition selection. This work focused on a tree structure with
composition on leaves, which is different from the compositional tree, where
composition exists at each node of the tree. Furthermore, it did not handle
boundary points (zeros or ones) in the data or cover asymptotic properties.

In this paper, we propose a regularized regression method to estimate the
association between a dependent variable and independent variables that have a
compositional tree structure. The regularization term is constructed from the
tree structure, which is assumed to be known, and designed to achieve sparsity
in both marginal and conditional effects from independent variables. Our model
is transformation-free and able to handle boundary points (zeros or ones) in
the data. We also establish consistency and model selection consistency of our
estimators building on results from Lee et al., (2015).

In the next section, we introduce an MRI data example. In Section 3, we define
the compositional tree and regression model. In Section 4, we present our
proposed method to estimate the regression coefficients. We evaluate the
performance of our proposed method through simulations in Section 5. The MRI
data application is provided in Section 6. Section 7 discusses future
directions.

##  2 Data example

Data used in the preparation of this article were obtained from the
Alzheimer’s Disease Neuroimaging Initiative (ADNI) database
(adni.loni.usc.edu). The ADNI was launched in 2003 as a public-private
partnership, led by Principal Investigator Michael W. Weiner, MD. The primary
goal of ADNI has been to test whether serial magnetic resonance imaging (MRI),
positron emission tomography (PET), other biological markers, and clinical and
neuropsychological assessment can be combined to measure the progression of
mild cognitive impairment (MCI) and early Alzheimer’s disease (AD). We focus
on the data set acquired by Liu et al., (2019) from the ADNI database.

The ADNI data set contains 819 subjects, which were diagnosed at the baseline
as cognitive normal (229 subjects), mild cognitive impairment (MCI, 402
subjects) or Alzheimer’s disease (AD, 188 subjects). For each subject, the
composite memory score, MRI data, and a wide variety of demographic,
behavioral and other non-imaging covariates were collected at several time
points. We focus on the composite memory score and MRI data. The composite
memory score was measured using data from the ADNI neuropsychological battery
and validated by Crane et al., (2012), with higher scores indicating better
memory. The MRI data consist of high-resolution T1-weighted images, which are
preprocessed and segmented through MRICloud (www.MRICloud.org, Mori et al.,,
2016), a public platform for multi-contrast imaging segmentation and
quantification. The preprocessing steps include orientation, inhomogeneity
correction, and histogram matching following with large deformation
diffeomorphic metric mapping (LDDMM). Post to preprocessing, segmentations
were obtained by fusing the multi-atlas labelling method (Tang et al.,, 2013).

For the MRI data, the five-level brain segmentation defines 321 brain regions,
which form a tree structure. At the first level of brain segmentation, the
whole brain is partitioned into 7 brain regions. At the second level of brain
segmentation, each of the 7 brain regions is further segmented into smaller
regions. At the finest level, there are 236 brain regions. Figure 1 displays
the tree structure of the 236 brain regions at the finest granularity. For
each brain region, we extracted its volume. Based on the five-level brain
segmentation procedure, the volume of a brain region is equal to the combined
volume of its subregions (after one segmentation). Furthermore, the combined
volume of brain regions at the finest level is equal to the ICV.

Structural MRI data have been commonly used to identify biomarkers of AD
(Vemuri and Jack,, 2010). For example, the density of neurofibrillary tangles
is an established pathological hallmark of AD, which can be reflected by MRI.
We hence focus on the association between memory declination, a common symptom
of AD, and brain volumes. For each subject, we use the MRI data acquired at
the initial screening, i.e. first time point, and the composite memory score
acquired on the same day as the MRI scan or the first post-imaging
measurement.

##  3 Model and Assumptions

###  3.1 Compositional tree

We first define tree structure using notation from graph theory. Let
$V=\\{X_{1},\dots,X_{p}\\}$ be a set of random variables with $0\leq X_{j}\leq
1$ for $j=1,\dots,p$. Let $E$ be a set of directed edges among
$X_{1},\dots,X_{p}$ with $E\subset\\{(X_{j}\rightarrow X_{k}):X_{j},X_{k}\in
V\\}$. For each edge $(X_{j}\rightarrow X_{k})\in E$, we call $X_{j}$ the
parent of $X_{k}$, and $X_{k}$ the child of $X_{j}$. $X_{j}$ is a leaf node if
it has no child and a root node if it has no parent. $X_{j}$ is an ancestor of
$X_{k}$ if the directed edges in $E$ can form a directed path from $X_{j}$ to
$X_{k}$, for example, $(X_{j}\rightarrow X_{s}),(X_{s}\rightarrow X_{k})\in
E$.

######  Definition 1.

$(V,E)$ forms a tree if (1) no $X_{j}$ is an ancestor of itself (i.e., $E$ not
containing any directed cycle), (2) $V$ contains only one root node and (3)
each $X_{j}$ has at most one parent.

In Definition 1, condition (1) defines a directed acyclic graph, and
conditions (2) and (3) are often made in defining a rooted tree in graph
theory. Figure 2 gives an example of tree structure with
$V=\\{X_{1},\dots,X_{10}\\}$ and $E=\\{(X_{10}\rightarrow
X_{9}),(X_{10}\rightarrow X_{8}),(X_{9}\rightarrow X_{1}),(X_{9}\rightarrow
X_{7}),(X_{7}\rightarrow X_{2}),(X_{7}\rightarrow X_{3}),(X_{8}\rightarrow
X_{4}),(X_{8}\rightarrow X_{5}),(X_{8}\rightarrow X_{6})\\}$.

$X_{10}$$X_{9}$$X_{8}$$X_{1}$$X_{7}$$X_{2}$$X_{3}$$X_{4}$$X_{5}$$X_{6}$ Figure
2: An example of a tree with $p=10$.

In our data example, we can define a tree given the hierarchical brain
segmentation. Let each $X_{j},j=1,\dots,321$, represent the volume of a brain
region $j$ and let $V$ be the set of all $X_{j}$. We regard $X_{j}$ as the
parent of $X_{k}$ if brain region $k$ is a subregion of $j$ defined by one-
step segmentation (i.e. there is no other subregion of $j$ that contains $k$).
If $X_{j}$ is the parent of $X_{k}$, we also call brain region $j$ the parent
of brain region $k$. Then the edge set $E$ is defined as the collection of all
parent-child relationships among brain regions and the (only) root node is the
ICV. For the $(V,E)$ defined above, condition (1) of Definition 1 holds by
construction, condition (2) follows because the root node is the ICV, and
condition (3) results from the fact that a region cannot be part of two
disjoint bigger regions.

Although we define the tree structure using notations of graph theory, we
emphasize that we do not associate the tree structure with conditional
independence or causal diagrams, such as in graphical probabilistic models
(Pearl,, 2009). Our tree solely represents the hierarchical structure among
$X_{1},\dots,X_{p}$ and is used to add compositional constraints, as described
below. Our goal is to study the association between an outcome of interest and
covariates $X_{j},j=1,\dots,p$, instead of the relationships among covariates.

Consider compositional constraints on $(X_{1},\dots,X_{p})$ complying with the
tree structure. Denoting $q$ as the number of leaf nodes, we can arrange the
indices of $X_{1},\dots,X_{p}$ such that the first $q$ variables
$(X_{1},\dots,X_{q})$ are the leaf nodes. For each $j=1,\dots,p$, let
$c(j)=\\{k:(X_{j}\rightarrow X_{k})\in E\\}$ denote the index set of children
of $X_{j}$ and let $|c(j)|$ denote the cardinality of $c(j)$ (i.e., the number
of children of $X_{j}$). We then have the following definition of a
compositional tree.

######  Definition 2.

Assume $(V,E)$ forms a tree and $X_{1},\dots,X_{q}$ are the leaf nodes. Then
$(V,E)$ forms a compositional tree if (1) $\sum_{j=1}^{q}X_{j}=1$ and (2)
$X_{j}=\sum_{k\in c(j)}X_{k}$ for each $j>q$.

In Definition 2, condition (1) imposes a compositional constraint on the leaf
nodes. Condition (2) requires that each parent node is equal to the summation
of its children. Conditions (1) and (2) together imply that the root node is a
constant 1. In the example shown in Figure 2, the constraints for a
compositional tree are $X_{7}=X_{2}+X_{3}$, $X_{8}=X_{4}+X_{5}+X_{6}$,
$X_{9}=X_{1}+X_{7}$ and $X_{10}=X_{8}+X_{9}=1$. For the case that $X_{j}$ has
only one child $X_{k}$, Definition 2 implies that $X_{j}=X_{k}$ and we hence
drop $X_{k}$ to avoid any replicate. In this paper, we assume the
compositional tree $(V,E)$ for a column vector of random variables
$\boldsymbol{X}=(X_{1},\dots,X_{p})^{t}$ is known.

Compositional trees generalize the structure on compositional data by allowing
more constraints on $\boldsymbol{X}$ and are less studied. Although
$X_{q+1},\dots,X_{p}$ are linear combinations of leaf nodes, they still
provide information on the structure of $\boldsymbol{X}$ and can help
interpret conditional effects (defined in Section 3.2 below). To simplify
notation, we say that $\boldsymbol{X}$ has a compositional tree structure if
the associated $(V,E)$ forms a compositional tree.

In our data example, brain regions $\boldsymbol{X}$ defined above have a
compositional tree structure. Since the summation of all leaf node volumes is
the ICV, then condition (1) of Definition 2 requires that the volumetric data
are normalized by the ICV such that each person has a total brain volume 1.
This is common practice in MRI analysis, since the ICV is typically only
meaningfully related to physical size. Alternative strategies remove
ventricular volumes and then study regional volumes relative to total brain
volume (i.e., studying tissue composition). We include the ventricular volumes
and normalize by ICV, since they are an important aspect of understanding
progressive tissue loss in a disorder like AD. Condition (2) of Definition 2
states that the volume of each brain region is equal to the combined volume of
all its children, which follows by definition that each brain region is
partitioned into subregions with no volume left undefined.

When dealing with compositional data, most current models work on a
transformed space, for example, Isometric logratio transformations (Egozcue et
al.,, 2003) or log ratio transformations (Papke and Wooldridge,, 1996).
Although such transformations provide convenience in estimation, they cannot
handle boundary values in $\boldsymbol{X}$ and add difficulty to
interpretation (Fiksel et al.,, 2020). We hence work on the original space
$\\{\boldsymbol{X}:X_{j}\geq
0,j=1,\dots,p;\sum_{j=1}^{p}X_{j}=1;X_{j}=\sum_{k\in
c(j)}X_{k},j=q+1,\dots,p\\}$.

For a compositional tree, the vector space spanned by $\boldsymbol{X}$ has
dimension (at most) $q<p$, which causes rank deficiency in many regression
models. An alternative way is to model $X_{j}=\sum_{k\in
c(j)}X_{k}+\varepsilon_{j}$, where $\varepsilon_{j}$ is an independent
Gaussian noise, following the method of Shojaie and Michailidis, (2010).
Although this method does not have the competing issue of rank deficiency, as
long as the covariance matrix of $(\varepsilon_{1},\dots,\varepsilon_{p})$ is
positive definite, we have $\varepsilon_{j}=0$ almost always, which violates
Gaussian modeling assumptions.

###  3.2 Linear model, parameter identifiability and interpretation

Let $Y$ be the outcome of interest. We assume the following linear model

$Y=\sum_{j=1}^{p}\beta_{j}X_{j}+\varepsilon=\boldsymbol{\beta}^{\top}\boldsymbol{X}+\varepsilon,$
(1)

where $\boldsymbol{\beta}=(\beta_{1},\dots,\beta_{p})^{\top}$ is a column
vector of unknown parameters, $\boldsymbol{X}$ has a compositional tree
structure with $q$ leaf nodes, and $\varepsilon\sim N(0,\sigma^{2})$ is
independent of $\boldsymbol{X}$. Since the root node of a compositional tree
is a constant 1 and included in $\boldsymbol{X}$, the intercept term is
omitted from model (1). For $i=1,\dots,n$, let
$(\boldsymbol{X}_{i},\varepsilon_{i})$ be independent, identically distributed
samples from the joint distribution of $(\boldsymbol{X},\varepsilon)$ and let
$Y_{i}=\boldsymbol{\beta}^{\top}\boldsymbol{X}_{i}+\varepsilon_{i}$.

Since $\boldsymbol{X}$ is rank deficient (with rank at most $q$),
$\boldsymbol{\beta}$ is not unique. Due to this fact, each
$\beta_{j},j=1,\dots,p$ is not interpretable without further assumptions. To
overcome this difficulty, we impose the following $p-q$ linear constraints on
$\boldsymbol{\beta}$:

$\sum_{k\in c(j)}\beta_{k}=0\textrm{ for all }j>q,$ (2)

which uniquely define a $\boldsymbol{\beta}\in\mathcal{S}$ (as shown in the
Supplementary Material). Linear constraints (2) require that, for each $X_{j}$
that is not a leaf node, the average effect of its children on $Y$ is $0$.
Then, each $\beta_{k}$ can be interpreted as the deviation effect of $X_{k}$
from the effect of its parent, $X_{j}$, on $Y$. To show this, consider the
following derivation using Definition 2:

$\beta_{j}X_{j}+\sum_{k\in
c(j)}\beta_{k}X_{k}=\beta_{j}=\left(\beta_{j}+\frac{1}{|c(j)|}\sum_{l\in
c(j)}\beta_{l}\right)X_{j}+\sum_{k\in
c(j)}\left(\beta_{k}-\frac{1}{|c(j)|}\sum_{l\in c(j)}\beta_{l}\right)X_{k},$

which implies that the average coefficient of children of $X_{j}$ can be
absorbed into the coefficient of $X_{j}$ and hence the remaining coefficients
of $X_{k},k\in c(j)$ are the deviations from $X_{j}$. By repeating this
procedure recursively from leaf nodes to the root node, we get all
coefficients satisfying linear constraints (2) with the desired
interpretation. For conciseness, $\beta_{k}$ is referred to as the
“conditional deviation effect” throughout, since its interpretation is
conditioning on the parent of $X_{k}$, i.e. the parent of $X_{k}$ held
constant.

Let $X_{p}$ denote the root node and $a(j)$ be the index set of ancestors of
$X_{j}$. Then the linear model (1) with constraints (2) can be formulated as:

$\displaystyle Y$
$\displaystyle=\beta_{p}+\sum_{j=1}^{q}\alpha_{j}X_{j}+\varepsilon=\beta_{p}+\boldsymbol{\alpha}^{\top}\boldsymbol{X}_{leaf}+\varepsilon,$
(3) $\displaystyle\quad\textrm{subject to}\quad\sum_{j=1}^{q}\alpha_{j}=0,$
(4)

where $\beta_{p}$ is the regression coefficient of the root node $X_{p}$ and
serves as the intercept,
$\boldsymbol{\alpha}=(\alpha_{1},\dots,\alpha_{q})^{\top}$ with
$\alpha_{j}=\beta_{j}+\sum_{k\in a(j)\setminus\\{p\\}}\beta_{k}$ and
$\boldsymbol{X}_{leaf}=(X_{1},\dots,X_{q})^{\top}$ is the vector of leaf
nodes. We assume that the only linear constraint on $\boldsymbol{X}_{leaf}$ is
$\sum_{j=1}^{q}X_{j}=1$, i.e., no component of $\boldsymbol{X}_{leaf}$ being a
linear combination of the others. The model (3) not only provides direct
interpretation of marginal associations between $Y$ and
$\boldsymbol{X}_{leaf}$ (which we introduce below), but is also useful for
estimating $\boldsymbol{\beta}$ in Section 4.

Compared to model (1), model (3) only uses the leaf nodes. Each $\alpha_{j}$
is the aggregation of conditional deviation effect of ancestors of $X_{j}$
excluding the root node. For each $j=1,\dots,q$, $\alpha_{j}$ can be
interpreted as the deviation effect of $X_{j}$ from the average effect of all
leaf nodes, referred to as the “marginal deviation effect” throughout for
conciseness. If $X_{j}$ is increased by $\delta$ at the expense of another
leaf node, $X_{k}$, i.e., $X_{k}$ decreased by $\delta$, then $Y$ is changed
by $\alpha_{j}-\alpha_{k}$. If $X_{j}$ is increased by $\delta$ at the expense
of all other leaf nodes evenly, i.e., $X_{k}$ decreased by $\delta/(q-1)$ for
all $k\leq q$ and $k\neq j$, then $Y$ is changed by
$\frac{q}{q-1}\delta\alpha_{j}$. Without the constraint (4), $\alpha_{j}$ is
not identifiable, since $\sum_{j=1}^{q}X_{j}=1$. However,
$\beta_{p}+\alpha_{j}$ would still be identifiable and a fact that we use for
estimation in Section 4. This statement is formally described in Proposition 1
below, which is proven in the Supplementary Material.

######  Proposition 1.

Assume that $\sum_{j=1}^{q}X_{j}=1$ and no component of
$\boldsymbol{X}_{leaf}$ is a linear combination of the others. If two sets of
parameters $(\beta_{p},\alpha_{1},\dots,\alpha_{q})$ and
$(\overline{\beta}_{p},\overline{\alpha}_{1},\dots,\overline{\alpha}_{q})$
both satisfy model (3), then
$\beta_{p}+\alpha_{j}=\overline{\beta}_{p}+\overline{\alpha}_{j}$ for each
$j=1,\dots,q$.

Compared to $\boldsymbol{\alpha}$, which is the marginal deviation effect,
$\boldsymbol{\beta}$ is the conditional deviation effect, offering flexibility
for interpreting various conditional effects. For example, in the tree
structure shown in Figure 2, $\beta_{2}+\beta_{7}$ represents the deviation
effect of $X_{2}$ from $X_{9}$, and $\beta_{7}+\beta_{9}$ represents the
marginal deviation effect of $X_{7}$ (i.e. conditioning on a constant
$X_{10}$).

In our data example, both $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ are
scientifically meaningful. The marginal deviation effect represents the effect
of the fractional volume of a leaf region on memory, while the conditional
deviation effect is the residual effect of the fractional volume of a brain
region on memory after removing the effect of its ancestors on memory.

##  4 Estimation

Let $\boldsymbol{\beta}^{*}$, $\boldsymbol{\alpha}^{*}$ denote the true
parameters that satisfy model (1) with constraints (2) and model (3) with
constraint (4) respectively. Our goal is to estimate $\boldsymbol{\alpha}^{*}$
and $\boldsymbol{\beta}^{*}$. Since $p$ and $q$ are potentially large ($p=321$
and $q=236$ for our data example), we propose a new regularization term based
on lasso for component selection (Lin et al.,, 2014) and fused lasso
(Tibshirani et al.,, 2005) and perform regularized regression to achieve
sparsity in both $\widehat{\boldsymbol{\alpha}}$ and
$\widehat{\boldsymbol{\beta}}$. In the method described below, we first
estimate $\boldsymbol{\alpha}^{*}$ using the generalized lasso Tibshirani and
Taylor, (2011) and then calculate $\widehat{\boldsymbol{\beta}}$ based on
$\widehat{\boldsymbol{\alpha}}$ by solving linear systems.

###  4.1 Regularization

For any $\boldsymbol{\beta}\in\mathbb{R}^{p}$ and
$\boldsymbol{\alpha}\in\mathbb{R}^{q}$, consider the regularization term

$P(\boldsymbol{\alpha},\boldsymbol{\beta},\eta)=\eta
P_{1}(\boldsymbol{\alpha})+(1-\eta)P_{2}(\boldsymbol{\beta}),$

where $\eta\in[0,1]$ is a tuning parameter adjusting the weight between two
$P_{1}(\boldsymbol{\alpha})$ and $P_{2}(\boldsymbol{\beta})$,

$\displaystyle P_{1}(\boldsymbol{\alpha})$
$\displaystyle=\sum_{j=1}^{q}\Big{|}\alpha_{j}-\frac{1}{q}\sum_{k=1}^{q}\alpha_{k}\Big{|},$
$\displaystyle P_{2}(\boldsymbol{\beta})$
$\displaystyle=\sum_{j=q+1}^{p}\sum_{s=1}^{|c(j)|-1}|\beta_{j_{s}}-\beta_{j_{s+1}}|,$

where $c(j)$ is the index set of children of $X_{j}$ with the elements in
$c(j)$ encoded as $j_{1},\dots,j_{|c(j)|}$. $P_{1}(\boldsymbol{\alpha})$
selects leaf nodes with non-zero marginal deviation effects. If
$\alpha_{j}=\frac{1}{q}\sum_{k=1}^{q}\alpha_{k}$, then changing $X_{j}$ at the
expense of all other leaf nodes evenly will not result in changes of $Y$. This
penalty is known as the lasso for component selection, which is also seen in
Wang and Zhao, (2017) for dealing with compositional data. In
$P_{2}(\boldsymbol{\beta})$, for each $X_{j}$ with $j>q$, we penalize the
difference among coefficients of its children using the fused lasso penalty.
If $|\beta_{j_{s}}-\beta_{j_{s+1}}|=0$ for all $s=1,\dots,|c(j)|-1$, which
means all children of $X_{j}$ have no conditional deviation effect, then the
component $\sum_{k\in c(j)}\beta_{k}X_{k}=\beta_{j_{1}}X_{j}$, resulting in a
sparse representation of linear model (1). Combined with the linear
constraints (2), the above case is also equivalent to $\beta_{k}=0$ for all
$k\in c(j)$. The following proposition gives some properties of
$P(\boldsymbol{\alpha},\boldsymbol{\beta},\eta)$.

######  Proposition 2.

Given linear constraints (2), there exists a matrix
$\boldsymbol{D}(\eta)\in\mathbb{R}^{(2q-1)\times q}$ such that
$P(\boldsymbol{\alpha},\boldsymbol{\beta},\eta)=||\boldsymbol{D}(\eta)\boldsymbol{\alpha}||_{1}$
and $\boldsymbol{D}(\eta)\boldsymbol{1}_{q}=\boldsymbol{0}_{q}$, where
$||\cdot||_{1}$ is the $L_{1}$-norm,
$\boldsymbol{1}_{q},\boldsymbol{0}_{q}\in\mathbb{R}^{q}$ are column vectors
with all entries 1, 0 respectively.

Proposition 2 implies that the penalty
$P(\boldsymbol{\alpha},\boldsymbol{\beta},\eta)$ can be formulated as a
function of $\boldsymbol{\alpha}$ and $\eta$, making it possible to perform
regularized regression based on model (3), which does not involve
$\boldsymbol{\beta}$. Furthermore, this penalty is invariant with respect to
constant change of $\boldsymbol{\alpha}$ (i.e.,
$||\boldsymbol{D}(\eta)\boldsymbol{\alpha}||_{1}=||\boldsymbol{D}(\eta)(\boldsymbol{\alpha}+C\boldsymbol{1}_{q})||_{1}$
for any $C\in\mathbb{R}$), which makes it equivalent to penalize on
$\boldsymbol{\alpha}+\beta_{q}\boldsymbol{1}_{q}$ as we do below. We prove
Proposition 2 and show how $\boldsymbol{D}(\eta)$ is constructed in the
Supplementary Material.

###  4.2 Estimating $\boldsymbol{\alpha}^{*}$

We estimate $\boldsymbol{\alpha}^{*}$ by
$\widehat{\boldsymbol{\alpha}}=\widehat{\widetilde{\boldsymbol{\alpha}}}-\boldsymbol{1}_{q}\boldsymbol{1}_{q}^{\top}\widehat{\widetilde{\boldsymbol{\alpha}}}$,
where

$\widehat{\widetilde{\boldsymbol{\alpha}}}=\arg\min_{\widetilde{\boldsymbol{\alpha}}}\frac{1}{n}\sum_{i=1}^{n}\left(Y_{i}-\widetilde{\boldsymbol{\alpha}}^{\top}\boldsymbol{X}_{leaf,i}\right)^{2}+\lambda||\boldsymbol{D}(\eta)\widetilde{\boldsymbol{\alpha}}||_{1}$
(5)

with $\boldsymbol{X}_{leaf,i}=(X_{i1},\dots,X_{iq})^{\top}$,
$\widetilde{\boldsymbol{\alpha}}=\boldsymbol{\alpha}+\beta_{p}\boldsymbol{1}_{q}$
and $\lambda>0$ being the tuning parameter. In equation (5),
$\widehat{\widetilde{\boldsymbol{\alpha}}}$ is an estimate of
$\widetilde{\boldsymbol{\alpha}}$, which is identifiable as discussed in
Section 3.2 and does not involve any linear constraints. Then,
$\widehat{\boldsymbol{\alpha}}$ is constructed by imposing the constraint (4),
i.e. centering $\widehat{\widetilde{\boldsymbol{\alpha}}}$. We note that the
regularization term
$\lambda||\boldsymbol{D}(\eta)\widetilde{\boldsymbol{\alpha}}||_{1}$ imposes
the desired sparsity on $\boldsymbol{\alpha}$, since
$\lambda||\boldsymbol{D}(\eta)\widetilde{\boldsymbol{\alpha}}||_{1}=\lambda||\boldsymbol{D}(\eta){\boldsymbol{\alpha}}||_{1}$
given Proposition 2.

For any $\boldsymbol{\alpha}\in\mathbb{R}^{q}$ and given $\eta\in[0,1]$, let
$\mathcal{S}(\boldsymbol{\alpha})$ be the support of
$\boldsymbol{D}(\eta)\boldsymbol{\alpha}$, i.e.,
$\mathcal{S}(\boldsymbol{\alpha})=\\{j\in\\{1,\dots,2q-1\\}:\boldsymbol{e}_{j}^{\top}\boldsymbol{D}(\eta)\boldsymbol{\alpha}\neq
0\\}$ with $\boldsymbol{e}_{j}\in\mathbb{R}^{2q-1}$ being a column vector with
the $j$-th entry 1 and the rest 0. Let
$\mathcal{M}=\\{\boldsymbol{\alpha}:\mathcal{S}(\boldsymbol{\alpha})\subset\mathcal{S}(\boldsymbol{\alpha}^{*})\\}$
denote the model subspace of interest. That is, for
$\boldsymbol{\alpha}\in\mathcal{M}$, an entry of
$\boldsymbol{D}(\eta)\boldsymbol{\alpha}$ is non-zero only if the
corresponding entry of $\boldsymbol{D}(\eta)\boldsymbol{\alpha}^{*}$ is non-
zero. The following theorem gives consistency and model selection consistency
of $\widehat{\boldsymbol{\alpha}}$, which is adapted from Corollary 4.2 of Lee
et al., (2015).

######  Theorem 1.

Given $\eta\in[0,1]$, we assume $\\{\boldsymbol{X}_{leaf,i}\\}_{i=1}^{n}$
satisfies restricted strong convexity (RSC) on $\mathcal{M}$ and
irrepresentability, which we define in the Supplementary Material. For
$\lambda=C_{1}\sigma\sqrt{\frac{\log q}{n}}$, $\widehat{\boldsymbol{\alpha}}$
is unique and, with probability at least $1-2/q$,

  1. 1.

(consistency)
$||\widehat{\boldsymbol{\alpha}}-\boldsymbol{\alpha}^{*}||_{2}\leq
C_{2}\sigma\sqrt{\frac{\log q}{n}}$,

  2. 2.

(model selection consistency) $\widehat{\boldsymbol{\alpha}}\in\mathcal{M}$,

where $||\cdot||_{2}$ is the $L_{2}$-norm and $C_{1},C_{2}$ are known
constants given in the Supplementary Material.

Theorem 1 implies that when $q$ and $n/\log(q)$ are large, then, with high
probability, our estimate $\widehat{\boldsymbol{\alpha}}$ is close to the
truth and does not contain false positives (non-zero effect of inactive
predictors with respect to $\boldsymbol{D}(\eta)$). The RSC assumption is
typically satisfied when $X_{leaf,i}$ follows a multivariate normal
distribution (Raskutti et al.,, 2010). The irrepresentability assumption
requires that the active predictors (with respect to $\boldsymbol{D}(\eta)$)
are not overly well-aligned with the inactive predictors. This is achieved
when the inactive predictors are orthogonal or nearly-orthogonal to the active
predictors (Lee et al.,, 2015). We provide a detailed description and
discussion of these assumptions in the Supplementary Material.

Given $\eta$, the optimization problem (5) can be solved by the genlasso
package (Tibshirani and Taylor,, 2011) in R software. To select the tuning
parameter $\lambda$, we propose to use the Akaike information criterion (AIC,
Akaike et al.,, 1998) or Bayesian information criterion (BIC, Schwarz,, 1978).
Let

$IC_{\gamma}(\eta,\lambda)=n\log\left\\{\sum_{i=1}^{n}\left(Y_{i}-\widehat{\widetilde{\boldsymbol{\alpha}}}^{\top}\boldsymbol{X}_{leaf,i}\right)^{2}\right\\}+\gamma\
\textrm{df}(\eta,\lambda),$

where $\gamma$ is a complexity factor, $\textrm{df}(\eta,\lambda)$ is the
effective number of parameters in $\widehat{\widetilde{\boldsymbol{\alpha}}}$.
$IC_{\gamma}(\eta,\lambda)$ refers to AIC if $\gamma=2$ and BIC if
$\gamma=\log(n)$. For any $\eta\in[0,1]$, define
$\widehat{\lambda}(\eta)=\arg\min_{\lambda\geq 0}IC_{\gamma}(\eta,\lambda)$.
We select the tuning parameters
$\widehat{\eta}=\arg\min_{\eta\in[0,1]}IC_{\gamma}(\eta,\widehat{\lambda}(\eta))$
and $\widehat{\lambda}=\widehat{\lambda}(\widehat{\eta})$. An alternative
method to tune parameters is cross-validation, but we do not consider it here,
since it would dramatically increase the computation complexity and performs
similarly to AIC.

###  4.3 Estimating $\boldsymbol{\beta}^{*}$

Given $\widehat{\boldsymbol{\alpha}}$, we calculate
$\widehat{\boldsymbol{\beta}}$ as follows. Since
$\alpha_{j}=\beta_{j}+\sum_{k\in a(j)}\beta_{k}-\beta_{p}$ for $j=1,\dots,q$,
we can construct a matrix $\boldsymbol{Q}_{1}\in\mathbb{R}^{q\times p}$ such
that $\boldsymbol{Q}_{1}\boldsymbol{\beta}=\boldsymbol{\alpha}$. Since
$\boldsymbol{\beta}$ also satisfies linear constraints (2), we can construct
another matrix $\boldsymbol{Q}_{2}\in\mathbb{R}^{(q-p)\times p}$ such that
$\boldsymbol{Q}_{2}\boldsymbol{\beta}=\boldsymbol{0}_{p-q}$. Denoting
$\boldsymbol{Q}=(\boldsymbol{Q}_{1}^{\top},\boldsymbol{Q}_{2}^{\top})^{\top}$,
then $\widehat{\boldsymbol{\beta}}$ is calculated by solving the linear system

$\boldsymbol{Q}\boldsymbol{\beta}=\left(\begin{array}[]{c}\widehat{\boldsymbol{\alpha}}\\\
\boldsymbol{0}_{p-q}\end{array}\right).$ (6)

The following theorem implies that $\widehat{\boldsymbol{\beta}}$ is uniquely
determined by $\widehat{\boldsymbol{\alpha}}$ (i.e. $\boldsymbol{Q}$ is
invertible) and is consistent and model selection consistent under the same
conditions as $\widehat{\boldsymbol{\alpha}}$.

######  Theorem 2.

Let $C_{1},C_{2},\lambda$ and $\mathcal{M}$ be the quantities defined in
Theorem 1. Given the same assumptions made in Theorem 1,
$\widehat{\boldsymbol{\beta}}$ is uniquely determined by
$\widehat{\boldsymbol{\alpha}}$ and, with probability at least $1-2/q$,

  1. 1.

(consistency)
$||\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{*}||_{2}\leq||\boldsymbol{Q}^{-1}||_{2}C_{2}\sigma\sqrt{\frac{\log
q}{n}}$,

  2. 2.

(model selection consistency)
$\boldsymbol{Q}_{1}\widehat{\boldsymbol{\beta}}\in\mathcal{M}$.

An alternative method to estimate $\boldsymbol{\beta}$ is solving a
constrained optimization problem following Lin et al., (2014):

$\displaystyle\widehat{\boldsymbol{\beta}}=$
$\displaystyle\arg\min_{\boldsymbol{\beta}}\frac{1}{n}\sum_{i=1}^{n}\left(Y_{i}-\boldsymbol{\beta}^{\top}\boldsymbol{X}_{i}\right)^{2}+\lambda
P(\alpha,\boldsymbol{\beta},\eta),$ subject to $\displaystyle\ \textrm{linear
constraints\leavevmode\nobreak\ (\ref{eq:beta-constraint})}.$

However, this method has to handle the rank deficiency of $\boldsymbol{X}$ and
$p-q$ linear constraints. If $q$ is much smaller than $p$, then the linear
constraints can be large, which may cause bias and increased computational
complexity. In our proposed method, these two issues are avoided by using two
steps to estimate $\boldsymbol{\beta}^{*}$ (first estimating
$\boldsymbol{\alpha}^{*}$ and then $\boldsymbol{\beta}^{*}$).

To the best of our knowledge, we are the first to study regularized regression
on a compositional tree and provide consistency and model selection
consistency. Kim and Xing, (2012) developed a tree-guided group lasso method,
but their goal was to analyze multi-response data and they did not consider
composition. Lin et al., (2014) used the lasso for component selection in
compositional data and their optimization problem is a special case of ours,
setting $\eta=1$. Wang and Zhao, (2017) proposed TASSO to perform penalized
regression on compositional data with a hierarchical tree structure, but they
did not consider compositional trees or provide asymptotic results. We compare
our method with lasso for component selection and TASSO in both simulations
and MRI data application below.

##  5 Simulation study

In this simulation study, we consider 4 data generating distributions, which
cover combinations of the following settings: a binary compositional tree or
the MRI-motivated compositional tree, and leaf or stem effects. A binary
compositional tree is a compositional tree where each parent has two children,
while the MRI-motivated compositional tree represents the same tree structure
as our data example (where a parent node may have more than two children).
Leaf effects stand for linear models where the true effects (non-zero
$\beta_{j}$) are only from nodes near the leaves, while stem effects mean that
true effects are only from nodes near the root. Different from the leaf
effects where both $\boldsymbol{\alpha}^{*}$ and $\boldsymbol{\beta}^{*}$ are
sparse, stem effects will lead to non-sparse $\boldsymbol{\alpha}^{*}$.

The first scenario (Scenario 1) has a binary compositional tree and leaf
effects. The tree structure is shown in Figure 3, where $p=255$, $q=128$ and
$\boldsymbol{X}_{leaf}=(X_{1},X_{2},\dots,X_{q})$. Letting $n=120$, we
independently generate $\boldsymbol{X}_{leaf,i},i=1,\dots,n$ by first
independently sampling $\widetilde{\boldsymbol{X}}_{leaf,i}$ from a
multivariate Gaussian distribution with mean $\boldsymbol{0}_{q}$ and variance
$\boldsymbol{\Sigma}=(\sigma_{ij})_{q\times q}$, where
$\sigma_{ij}=0.2^{|i-j|}$ is the $i$-th row $j$-th column entry of
$\boldsymbol{\Sigma}$, and then defining
$\boldsymbol{X}_{leaf,i}=\widetilde{\boldsymbol{X}}_{leaf,i}/\boldsymbol{1}_{q}^{\top}\widetilde{\boldsymbol{X}}_{leaf,i}$
to satisfy the composition condition. For $j>q$, we generate $X_{ij}$
following the definition of compositional tree using
$\boldsymbol{X}_{leaf,i}$. We define, for $i=1,\dots,n$,

$Y_{i}=3+X_{i,1}-X_{i,2}+X_{i,129}-X_{i,130}+\varepsilon=3+2X_{i,1}-X_{i,3}-X_{i,4}+\varepsilon,$

where $\varepsilon_{i}$ is an independent sample from $N(0,\sigma^{2})$ and
$\sigma^{2}$ is chosen such that
$Var(\boldsymbol{X}^{\top}\boldsymbol{\beta})=Var(\varepsilon)$. This model
only involves the left bottom corner in the tree shown in Figure 3. The non-
zero conditional deviation effects are
$\beta_{1}^{*}=\beta_{129}^{*}=1,\beta_{2}^{*}=\beta_{130}^{*}=-1$ and the
non-zero marginal deviation effects are
$\alpha_{1}^{*}=2,\alpha_{3}^{*}=\alpha_{4}^{*}=-1$.

The second scenario (Scenario 2) has a binary compositional tree and stem
effects, where the binary compositional tree and
$\boldsymbol{X}_{i},i=1,\dots,n$ is the same as in Scenario 1. For the stem
effect, we define

$Y_{i}=3+X_{i,249}-X_{i,250}+X_{i,253}-X_{i,254}+\varepsilon=3+2\sum_{j=1}^{32}X_{ij}-\sum_{j=65}^{128}X_{ij}+\varepsilon_{i},$

where $\varepsilon_{i}$ is defined in the same way as in Scenario 1. Unlike
Scenario 1, this model only involves the conditional deviation effects from
the top part in the tree (nodes near the root), which are
$\beta_{249}^{*}=\beta_{253}^{*}=1,\beta_{250}^{*}=\beta_{254}^{*}=-1$.
Furthermore, the marginal deviation effect $\boldsymbol{\alpha}$ is no longer
sparse because $\alpha_{j}^{*}=2$ for $j=1,\dots,32$ and $\alpha_{j}^{*}=-1$
for $j=65,\dots,128$.

$X_{255}$$X_{253}$$X_{254}$$X_{249}$$X_{250}$$X_{251}$$X_{252}$$\cdots$$\cdots$$\cdots$$\cdots$$\cdots$$\cdots$$\cdots$$\cdots$$X_{129}$$X_{130}$$X_{1}$$X_{2}$$X_{3}$$X_{4}$$\cdots$$\cdots$
Figure 3: The binary compositional tree considered in Scenarios 1 and 2 of the
simulation study with $p=255$ and $q=128$.

In the third scenario (Scenario 3), we consider the MRI-motivated
compositional tree and leaf effects. The MRI-motivated compositional tree is
shown in Figure 1, where $p=321$ and $q=236$. For our MRI data example,
$n=819$ and we denote the empirical distribution of
$(\widetilde{\boldsymbol{X}}_{1},\dots,\widetilde{\boldsymbol{X}}_{n})$ by
$F_{n}$, where $\widetilde{\boldsymbol{X}}_{i}$ contains the fractional brain
volumetric data of participant $i$. Let $\boldsymbol{X}_{i},i=1,\dots,n$ be
independent samples from $F_{n}$. We model, for $i=1,\dots,n$,

$Y_{i}=3+3\ X_{i,\ \textrm{SFG-L}}-2\ X_{i,\ \textrm{SFG-PFC-L}}-X_{i,\
\textrm{SFG-pole-L}}+\varepsilon,$

where SFG-L, SFG-PFC-L and SFG-pole-L are all leaf nodes and subregions of the
superior frontal gyrus left hemisphere and $\varepsilon$ is as defined in
Scenario 1. In this model, we have
$\beta_{\textrm{SFG-L}}^{*}=\alpha_{\textrm{SFG-L}}^{*}=3$,
$\beta_{\textrm{SFG-PFC-L}}^{*}=\alpha_{\textrm{SFG-PFC-L}}^{*}=-2$ and
$\beta_{\textrm{SFG-pole-L}}^{*}=\alpha_{\textrm{SFG-pole-L}}^{*}=-1$.

In the last scenario (Scenario 4), we consider the MRI-motivated compositional
tree again but with stem effects. We use the same compositional tree and
$\boldsymbol{X}_{i},i=1,\dots,n$ as in Scenario 3. Let

$Y=3+X_{i,\ \textrm{Telencephalon-L}}-X_{i,\
\textrm{Telencephalon-R}}+\varepsilon,$

where Telencephalon-L and Telencephalon-R represent telencephalon located in
the left and right hemisphere respectively and both are children of the ICV.
Different from Scenario 3, this model has
$\beta_{\textrm{Telencephalon-L}}^{*}=-1$,
$\beta_{\textrm{Telencephalon-R}}^{*}=1$ and 200 non-zero entries in
$\boldsymbol{\alpha}^{*}$.

For each of the 4 scenarios, we simulate $m=1000$ data sets and compare our
proposed method with lasso for component selection (abbreviated as CLASSO
throughout, Lin et al.,, 2014) and TASSO (Wang and Zhao,, 2017). For TASSO, we
use their default settings to estimate $\boldsymbol{\alpha}^{*}$ and calculate
$\widehat{\boldsymbol{\beta}}_{TASSO}$ by solving equation (6). The only
difference is that the natural log-transformation is not performed, as
described in Section 3. Since CLASSO is a special case of our proposed method,
we calculate $\widehat{\boldsymbol{\beta}}_{CLASSO}$ following the same
procedure as our method setting $\eta=1$. For all three methods, we use AIC or
BIC to select the tuning parameters. The following metrics are used to compare
their performances: (1) sensitivity, defined as $|\\{j:\widehat{\beta}_{j}\neq
0,{\beta}_{j}^{*}\neq 0\\}|/|\\{j:{\beta}_{j}^{*}\neq 0\\}|$, (2) specificity,
defined as
$|\\{j:\widehat{\beta}_{j}={\beta}_{j}^{*}=0\\}|/|\\{j:{\beta}_{j}^{*}=0\\}|$
and (3) sum squared error (SSE), defined as
$||\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{*}||_{2}^{2}$. For each of
the above metrics, we report its average and standard deviation over the $m$
data sets. Since Scenarios 1 and 2 have $n<q$, a small $L_{2}$-penalty
(0.0001) is added when solving the optimization problem (5).

Table 1 gives the simulation results for Scenarios 1-4. In Scenarios 1 and 3,
the true parameter $\boldsymbol{\alpha}^{*}$ is sparse and all three methods
perform well, as expected. For our proposed method, the tuning parameter,
$\eta$, is near 0.5, indicating regularization terms on both
$\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ help penalize. The
regularization term on $\boldsymbol{\beta}$ is not shown in TASSO and CLASSO
and hence leads to the slightly better performance of our method. In Scenarios
2 and 4, since the true parameter, $\boldsymbol{\alpha}^{*}$, is not sparse,
our proposed method outperforms the other two methods on all performance
metrics. In such cases, the $L_{1}$ penalty on $\boldsymbol{\alpha}$ does not
help. Hence, TASSO and CLASSO tend to over-penalize (low sensitivity, high
specificity) or under-penalize (high sensitivity, low specificity) on
$\boldsymbol{\alpha}$, either leading to high SSEs. In contrast, our proposed
method always selects $\eta=0$ under BIC tuning, implying it only penalizes on
differences of conditional deviation effects. Across all 4 scenarios, our
method has high accuracy. Compared to AIC, BIC tends to perform as well or
better mirroring the simulation results of Wang and Zhao, (2017). Hence BIC
tuning is used in the MRI data application. In Scenarios 1 and 2, since $n$ is
smaller than $q$ and $p$, all three methods have larger SSEs compared to
Scenarios 3 and 4. Our proposed method with BIC tuning, however, remains
accurate, suggesting its ability to deal with high-dimensional data.

In the Supplementary Material, we provide an additional simulation study,
where noisier data is used to stress-test the method. In particular,
$\sigma^{2}$ is set such that
$Var(\varepsilon)=10Var(\boldsymbol{X}^{\top}\boldsymbol{\beta})$ while all of
the other settings of Scenarios 1-4 are kept constant. In this simulation, all
methods perform worse, because of the weaker signal relative to the noise. Our
method, however, still outperforms TASSO and CLASSO and the findings described
in other simulations still hold. In addition, simulation results for Scenarios
1 and 2 setting $n=1000$ are provided, which show similarly good relative
performance.

Table 1: Simulation results for Scenarios 1-4 comparing our method, TASSO and
CLASSO.

| Method | Tuning | Sensitivity | Specificity | SSE | $\eta$
---|---|---|---|---|---|---
| Our method | AIC | 1(0.02) | 0.04(0.02) | 20.45(4.94) | 0.31(0.37)
|  | BIC | 0.97(0.12) | 0.96(0.08) | 0.8(1.13) | 0.49(0.17)
Scenario | TASSO | AIC | 1(0.02) | 0.05(0.02) | 18.64(4.58) | -
1 |  | BIC | 0.96(0.18) | 0.96(0.08) | 1.21(1.4) | -
| CLASSO | AIC | 1(0) | 0.04(0.02) | 19.41(4.81) | -
|  | BIC | 0.98(0.14) | 0.91(0.11) | 1.22(1.62) | -
| Our method | AIC | 1(0) | 0.01(0.01) | 881.89(205.73) | 0.37(0.38)
|  | BIC | 1(0) | 0.99(0.06) | 2(27.58) | 0(0.03)
Sceinaro | TASSO | AIC | 1(0) | 0.01(0.01) | 846.65(197.18) | -
2 |  | BIC | 0.04(0.19) | 0.96(0.19) | 28.96(127.08) | -
| CLASSO | AIC | 1(0) | 0.01(0.01) | 871.57(204.35) | -
|  | BIC | 0.42(0.48) | 0.89(0.26) | 57.89(181.14) | -
| Our method | AIC | 1(0) | 0.85(0.1) | 0.77(0.69) | 0.45(0.33)
|  | BIC | 1(0) | 0.96(0.03) | 0.78(0.57) | 0.4(0.24)
Scenario | TASSO | AIC | 1(0) | 0.92(0.05) | 1.68(0.44) | -
3 |  | BIC | 1(0) | 0.97(0.04) | 1.58(0.77) | -
| CLASSO | AIC | 1(0) | 0.79(0.07) | 2.41(0.4) | -
|  | BIC | 1(0) | 0.9(0.04) | 2.49(0.4) | -
| Our method | AIC | 1(0) | 0.93(0.09) | 0.9(3.32) | 0(0.01)
|  | BIC | 1(0) | 0.98(0.01) | 0.28(0.37) | 0(0)
Scenario | TASSO | AIC | 0.92(0.2) | 0.35(0.11) | 34.92(22.27) | -
4 |  | BIC | 0.48(0.46) | 0.84(0.12) | 13.12(10.08) | -
| CLASSO | AIC | 1(0.02) | 0.25(0.05) | 22.3(3.07) | -
|  | BIC | 0.99(0.11) | 0.59(0.12) | 14.17(4.72) | -

##  6 MRI data application

We applied our proposed method to the data example introduced in Section 2.
The outcome $Y$ is the composite memory score while $\boldsymbol{X}$ is the
brain volumes resulting from the five-level brain segmentation. Since the
simulation study shows that BIC outperforms AIC on the compositional tree of
the data example, we used BIC to tune the hyperparameters and obtained
$\widehat{\eta}=0.405$.

Our method identified 77 non-zero marginal deviation effects
($\widehat{\boldsymbol{\alpha}}$) from the 236 leaf brain regions. Because of
the composition property, each $\widehat{\alpha}_{j}$ can be other brain
regions as described in Section 3. Table 2 displays the 10 largest effects,
which accounts for 48% of $||\widehat{\boldsymbol{\alpha}}||_{1}$. Among the
10 largest effects, Hippo-L represents hippocampus in the left hemisphere,
which is a limbic subregion and whose atrophy is well established and studied
in the progression of AD (Pini et al.,, 2016). InferiorLV-L is the inferior
pars of the left lateral ventricle (LV). Evidence has shown that its
enlargement is related to MCI and AD (Nestor et al.,, 2008). The preponderance
at CSF structures highlights progressive atrophy in AD. Amyg-R stands for
amygdala in the right hemisphere. A recent study (Poulin et al.,, 2011) on
this region suggested that “the magnitude of amygdala atrophy is comparable to
that of the hippocampus in the earliest clinical stages of AD, and is related
to global illness severity.” SOG-L, MOG-L and IOG-L represent the left
superior, middle and inferior occipital gyri, respectively, and are identified
as a group (same marginal deviation effect) by our method. Cu-L and LG-L are
cuneus and lingual gyrus in the left occipital region respectively and also
have the same marginal deviation effect. Although the occipital subregions
have opposite signs of marginal deviation effects, their conditional deviation
effects (i.e. $\beta_{j}$) cancel off when combined, resulting in a positive
marginal deviation effect (25.06) of the occipital region on memory. Holroyd
et al., (2000) showed that occipital atrophy is associated with visual
hallucinations (the most common type of hallucination) in AD. However, less is
known about the different roles of occipital subregions in AD. SylParieSul-L
represents sylvian parietal sulcus in the left hemisphere. To the best of our
knowledge, its enlargement is associated with progression of AD (Liu et al.,,
2012), which is contrary to our finding. However, we note that this region is
also identified with positive marginal deviation effect by TASSO (43.87),
which may suggest a false-positive result of the variable selection methods or
a special structure of the data set. MTG-L stands for the left middle temporal
gyrus. Its atrophy has been associated with AD (Pini et al.,, 2016). However,
it is important to emphasize that these results are exploratory in nature,
since the method investigates a large possible collection of potential
relationships and we did not pre-register any specific hypotheses.

Table 2: Top 10 regression coefficients ($\alpha_{j}$) of MRI application. ROI | $\alpha_{j}$ | ROI | $\alpha_{j}$
---|---|---|---
Hippo-L | 518.64 | IOG-L | 90.32
InferiorLV-R | -261.88 | Cu-L | -70.34
Amyg-R | 172.73 | LG-L | -70.34
SOG-L | 90.32 | SylParieSul-L | 69.06
MOG-L | 90.32 | MTG-L | 62.14

The 77 marginal deviation effects are aggregations of 109 conditional
deviation effects ($\beta$). For the 10 largest marginal deviation effects, we
decomposed them into conditional deviation effects using the definition of
$\boldsymbol{\alpha}$ (Section 3) and displayed the results in Figure 4. All
10 effects are from CSF and telencephalon. The effects from the ventricle are
negative and the effects form the limbic region are positive, both of which
are consistent with existing scientific findings (Pini et al.,, 2016; Nestor
et al.,, 2008). Complete results for marginal and conditional deviation
effects are given in the Supplementary Material.

In addition to our proposed method, we also ran TASSO and CLASSO with BIC
tuning. Twenty non-zero marginal deviation effects and 73 non-zero conditional
deviation effects were identified by all three methods, which include brain
regions in the ventricles and temporal lobe, although the magnitude of these
effects differs substantially among methods. Especially, the Amyg-R (right
hemisphere amygdala) region is only identified by our method. Compared with
our method, TASSO and CLASSO identify fewer marginal effects (40 non-zero
entries in $\widehat{\boldsymbol{\alpha}}_{TASSO}$ and 27 non-zero entries in
$\widehat{\boldsymbol{\alpha}}_{CLASSO}$), but they have larger BIC (5114 for
our proposed method, 5134 for TASSO and 5120 for CLASSO), indicating a larger
residual error. In addition, our proposed method tends to group effects
together, e.g. the left hemisphere occipitial subregions, which can facilitate
the interpretation of marginal and conditional deviation effects. For all
three methods, the effects from left and right hemispheres are generally not
equal, potentially suggesting a laterally asymmetric correlation between
volume and memory.

Whole BrainTelencephalon-LTelencephalon-RCerebralCortex-LCerebralNucli-RAmyg-
RTemporal-LLimbic-LHippo-LOccipital-LMTG-L SOG-L MOG-L IOG-L Cu-L LG-L
CSFVentricleSulcus-LSylvianFissureExt-LLV-RInferiorLV-RSylParieSul-L(a)(b)

Figure 4: Conditional deviation effects $\boldsymbol{\beta}$ related to the 10
largest marginal deviation effects $\boldsymbol{\alpha}$. Suffix “-L” or “-R”
refers to the left or right hemisphere of the brain, respectively. Red (blue)
color represents positive (negative) sign of $\beta$ with darker color
indicating a larger value of $|\beta|$. Panels (a) and (b) show the aggregated
conditional effects in CSF and telencephalon in 3-dimensional template brain
space.

##  7 Discussion

The linear model in Section 3 also allows for including additional covariates,
in addition to covariates associated with the compositional tree. However,
when interaction terms are added, the linear constraints (2) can only handle
interactions between additional covariates and the whole compositional tree.

In our method, for estimating $\boldsymbol{\alpha}^{*}$, we assume that no
components of the leaf nodes are linear combinations of the others such that
$\widetilde{\boldsymbol{\alpha}}^{*}$ is identifiable. This assumption
generally holds if no further linear constraints are made on the leaf nodes.
When this assumption is not true, one can add a small $L_{2}$-penalty to the
right side of equation (5) and run the model, otherwise unmodified. In this
case, point estimates of $\boldsymbol{\alpha}^{*}$ and
$\boldsymbol{\beta}^{*}$ may be biased because of the $L_{2}$-penalty.

Our proposed method also assumes that the outcome is continuous. If the
outcome is binary or a count, then relatively minor modifications could use
generalized linear models. However, since the loss function is no longer
linear, how to consistently estimate $\boldsymbol{\alpha}^{*}$ with
generalized lasso penalty remains future research.

## Acknowledgements

Data collection and sharing for this project was funded by the Alzheimer’s
Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant
U01 AG024904) and DOD ADNI (Department of Defense award number
W81XWH-12-2-0012). ADNI is funded by the National Institute on Aging, the
National Institute of Biomedical Imaging and Bioengineering, and through
generous contributions from the following: AbbVie, Alzheimer’s Association;
Alzheimer’s Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.;
Biogen; Bristol-Myers Squibb Company; CereSpir, Inc.; Cogstate; Eisai Inc.;
Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La
Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE
Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Research &
Development, LLC.; Johnson & Johnson Pharmaceutical Research & Development
LLC.; Lumosity; Lundbeck; Merck & Co., Inc.; Meso Scale Diagnostics, LLC.;
NeuroRx Research; Neurotrack Technologies; Novartis Pharmaceuticals
Corporation; Pfizer Inc.; Piramal Imaging; Servier; Takeda Pharmaceutical
Company; and Transition Therapeutics. The Canadian Institutes of Health
Research is providing funds to support ADNI clinical sites in Canada. Private
sector contributions are facilitated by the Foundation for the National
Institutes of Health (www.fnih.org). The grantee organization is the Northern
California Institute for Research and Education, and the study is coordinated
by the Alzheimer’s Therapeutic Research Institute at the University of
Southern California. ADNI data are disseminated by the Laboratory for Neuro
Imaging at the University of Southern California.

## References

  * Akaike et al., (1998) Akaike, H., Parzen, E., Tanabe, K., and Kitagawa, G. (1998).  Selected papers of hirotugu akaike.  Springer Science & Business Media.
  * Crane et al., (2012) Crane, P. K., Carle, A., Gibbons, L. E., Insel, P., Mackin, R. S., Gross, A., Jones, R. N., Mukherjee, S., Curtis, S. M., Harvey, D., Weiner, M., Mungas, D., and for the Alzheimer’s Disease Neuroimaging Initiative (2012).  Development and assessment of a composite score for memory in the Alzheimer’s Disease Neuroimaging Initiative (ADNI).  Brain imaging and behavior, 6(4):502–516.
  * Egozcue et al., (2003) Egozcue, J. J., Pawlowsky-Glahn, V., Mateu-Figueras, G., and Barcelo-Vidal, C. (2003).  Isometric logratio transformations for compositional data analysis.  Mathematical Geology, 35(3):279–300.
  * Fiksel et al., (2020) Fiksel, J., Zeger, S., and Datta, A. (2020).  A transformation-free linear regression for compositional outcomes and predictors.  arXiv preprint arXiv:2004.07881.
  * Holroyd et al., (2000) Holroyd, S., Shepherd, M. L., and Downs III, J. H. (2000).  Occipital atrophy is associated with visual hallucinations in Alzheimer’s disease.  The Journal of neuropsychiatry and clinical neurosciences, 12(1):25–28.
  * Kim and Xing, (2012) Kim, S. and Xing, E. P. (2012).  Tree-guided group lasso for multi-response regression with structured sparsity, with an application to eQTL mapping.  The Annals of Applied Statistics, 6(3):1095–1117.
  * Lee et al., (2015) Lee, J. D., Sun, Y., and Taylor, J. E. (2015).  On model selection consistency of regularized m-estimators.  Electronic Journal of Statistics, 9(1):608–642.
  * Leite, (2016) Leite, M. L. C. (2016).  Applying compositional data methodology to nutritional epidemiology.  Statistical methods in medical research, 25(6):3057–3065.
  * Lin et al., (2014) Lin, W., Shi, P., Feng, R., and Li, H. (2014).  Variable selection in regression with compositional covariates.  Biometrika, 101(4):785–797.
  * Liu et al., (2019) Liu, C.-F., Padhy, S., Ramachandran, S., Wang, V. X., Efimov, A., Bernal, A., Shi, L., Vaillant, M., Ratnanather, J. T., and Faria, A. V. (2019).  Using deep siamese neural networks for detection of brain asymmetries associated with Alzheimer’s disease and mild cognitive impairment.  Magnetic resonance imaging, 64:190–199.
  * Liu et al., (2012) Liu, T., Lipnicki, D. M., Zhu, W., Tao, D., Zhang, C., Cui, Y., Jin, J. S., Sachdev, P. S., and Wen, W. (2012).  Cortical gyrification and sulcal spans in early stage Alzheimer’s disease.  PloS one, 7(2):e31083.
  * Ma and Zhang, (2020) Ma, X. and Zhang, P. (2020).  Quantile regression for compositional covariates.  arXiv preprint arXiv:2006.00789.
  * Mori et al., (2016) Mori, S., Wu, D., Ceritoglu, C., Li, Y., Kolasny, A., Vaillant, M. A., Faria, A. V., Oishi, K., and Miller, M. I. (2016).  Mricloud: delivering high-throughput MRI neuroinformatics as cloud-based software as a service.  Computing in Science & Engineering, 18(5):21–35.
  * MSCI, (2020) MSCI (2020).  Global industry classification standard (gics) methodology.  https://www.msci.com/documents/1296102/11185224/GICS+Methodology+2020.pdf/.
  * Mullahy, (2015) Mullahy, J. (2015).  Multivariate fractional regression estimation of econometric share models.  Journal of Econometric Methods, 4(1):71–100.
  * Nestor et al., (2008) Nestor, S. M., Rupsingh, R., Borrie, M., Smith, M., Accomazzi, V., Wells, J. L., Fogarty, J., Bartha, R., and Alzheimer’s Disease Neuroimaging Initiative (2008).  Ventricular enlargement as a possible measure of Alzheimer’s disease progression validated using the Alzheimer’s Disease Neuroimaging Initiative database.  Brain, 131(9):2443–2454.
  * Papke and Wooldridge, (1996) Papke, L. E. and Wooldridge, J. M. (1996).  Econometric methods for fractional response variables with an application to 401 (k) plan participation rates.  Journal of applied econometrics, 11(6):619–632.
  * Pawlowsky-Glahn and Egozcue, (2006) Pawlowsky-Glahn, V. and Egozcue, J. J. (2006).  Compositional data and their analysis: an introduction.  Geological Society, London, Special Publications, 264(1):1–10.
  * Pearl, (2009) Pearl, J. (2009).  Causality.  Cambridge university press.
  * Pini et al., (2016) Pini, L., Pievani, M., Bocchetta, M., Altomare, D., Bosco, P., Cavedo, E., Galluzzi, S., Marizzoni, M., and Frisoni, G. B. (2016).  Brain atrophy in Alzheimer’s disease and aging.  Ageing research reviews, 30:25–48.
  * Poulin et al., (2011) Poulin, S. P., Dautoff, R., Morris, J. C., Barrett, L. F., Dickerson, B. C., and Alzheimer’s Disease Neuroimaging Initiative (2011).  Amygdala atrophy is prominent in early Alzheimer’s disease and relates to symptom severity.  Psychiatry Research: Neuroimaging, 194(1):7–13.
  * Raskutti et al., (2010) Raskutti, G., Wainwright, M. J., and Yu, B. (2010).  Restricted eigenvalue properties for correlated Gaussian designs.  The Journal of Machine Learning Research, 11:2241–2259.
  * Schwarz, (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  The annals of statistics, 6(2):461–464.
  * Shojaie and Michailidis, (2010) Shojaie, A. and Michailidis, G. (2010).  Penalized likelihood methods for estimation of sparse high-dimensional directed acyclic graphs.  Biometrika, 97(3):519–538.
  * Tang et al., (2013) Tang, X., Oishi, K., Faria, A. V., Hillis, A. E., Albert, M. S., Mori, S., and Miller, M. I. (2013).  Bayesian parameter estimation and segmentation in the multi-atlas random orbit model.  PloS one, 8(6):e65591.
  * Tibshirani et al., (2005) Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., and Knight, K. (2005).  Sparsity and smoothness via the fused lasso.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(1):91–108.
  * Tibshirani and Taylor, (2011) Tibshirani, R. J. and Taylor, J. (2011).  The solution path of the generalized lasso.  The Annals of Statistics, 39(3):1335–1371.
  * Vemuri and Jack, (2010) Vemuri, P. and Jack, C. R. (2010).  Role of structural MRI in Alzheimer’s disease.  Alzheimer’s research & therapy, 2(4):23.
  * Wang and Zhao, (2017) Wang, T. and Zhao, H. (2017).  Structured subcomposition selection in regression and its application to microbiome data analysis.  The Annals of Applied Statistics, 11(2):771–791.