# Correlations in the Electronic Structure of van der Waals NiPS3 Crystals: An
X-Ray Absorption and Resonant Photoelectron Spectroscopy Study

Mouhui Yan  [  Yichen Jin  [  Zhicheng Wu  [  Arshak Tsaturyan  [  Anna
Makarova  [  Dmitry Smirnov  [  Elena Voloshina  [<EMAIL_ADDRESS>Yuriy Dedkov  [<EMAIL_ADDRESS>
###### Abstract

The electronic structure of high-quality van der Waals NiPS3 crystals was
studied using near-edge x-ray absorption spectroscopy (NEXAFS) and resonant
photoelectron spectroscopy (ResPES) in combination with density functional
theory (DFT) approach. The experimental spectroscopic methods, being element
specific, allow to discriminate between atomic contributions in the valence
and conduction band density of states and give direct comparison with the
results of DFT calculations. Analysis of the NEXAFS and ResPES data allows to
identify the NiPS3 material as a charge-transfer insulator. Obtained
spectroscopic and theoretical data are very important for the consideration of
possible correlated-electron phenomena in such transition-metal layered
materials, where the interplay between different degrees of freedom for
electrons defines their electronic properties, allowing to understand their
optical and transport properties and to propose further possible applications
in electronics, spintronics and catalysis.

SHU] Department of Physics, Shanghai University, 200444 Shanghai, China SHU]
Department of Physics, Shanghai University, 200444 Shanghai, China SHU]
Department of Physics, Shanghai University, 200444 Shanghai, China SFedU]
Institute of Physical and Organic Chemistry, Southern Federal University,
344090 Rostov on Don, Russia FUB] Institut für Chemie und Biochemie, Freie
Universität Berlin, Arnimallee 22,
14195 Berlin, Germany TUDD] Institut für Festkörper- und Materialphysik,
Technische Universität Dresden,
01069 Dresden, Germany SHU] Department of Physics, Shanghai University, 200444
Shanghai, China [FUB] Institut für Chemie und Biochemie, Freie Universität
Berlin, Arnimallee 22,
14195 Berlin, Germany SHU] Department of Physics, Shanghai University, 200444
Shanghai, China [FUB] Institut für Chemie und Biochemie, Freie Universität
Berlin, Arnimallee 22,
14195 Berlin, Germany

Layered materials, like graphene, silicene or transition metal dichalcogenides
(TMDs) and different heterosystems on their basis have attracted much
attention in the last years 1, 2, 3, 4, 5, 6. These studies are motivated by
the hope on the possible new applications of these materials in different
areas, where their unique electronic and transport properties might bring new
functionalities. These layered materials can be relatively easily produced in
the high quality form either as bulk materials or as a monolayer (or
multilayer) on different substrates and then isolated for further studies or
applications 7, 8, 9, 10, 11, 12. Such high-quality layers and heterosystems
made from layered 2D materials allow to perform studies of different
intriguing collective quantum phenomena, which appear in low-dimensional
systems, like, for example, superconductivity for twisted graphene bi-layers
13, 14, charge density waves in TMDs 15, 16, and coherent excitonic
excitations in antiferromagnetic transition metal trichalcogenide (TMT) NiPS3
17. For example, recently, different exciting phenomena in 2D materials were
extended to the observation of the long-range Ising-type antiferromagnetic
(AFM) order in the TMT FePS3 monolayer with a Néel temperature ($T_{N}$) of
$118$ K 18, 19. Following these works, a layer-dependent ferromagnetic (FM)
order in CrI3 with a Curie temperature ($T_{C}$) of $45$ K was reported 20 and
then the FM magnetic order has been observed for single atomic layers of
Cr2Ge2Te6 ($T_{C}\sim 40$ K) 21, Fe3GeTe2 ($T_{C}\sim 220$ K) 22, and MnSe2
23. These inspiring discoveries of magnetic order in 2D monolayers open up
exciting opportunities for the application of these properties in new
transport and optical devices and for the exploration of intrinsic 2D Ising-
type (anti)ferromagnetism.

The above mentioned layered transition metal trichalcogenides (TMPX3; TM: $3d$
TMs, like Fe, Co, or Ni; X: chalcogen, like S or Se) came to the focus of the
studies of 2D materials only recently 24, 25, 26. From the crystallographic
point of view they can be seen as layers of TMS2 where one third of TM atoms
is replaced by P–P dimers perpendicular to the layer (Fig. 1(a)) and in the 3D
structure they crystalize into the $C2/m$ or $R\bar{3}$ space groups 27, 28.
Most of these compounds are AFM wide-gap semiconductors 29, 30, 28 and they
were proposed for different applications like gas sensors, optical sensors,
effective materials in water splitting applications and in (opto)spintronics
31, 32, 26. Surprisingly, at the present time, the main focus in the studies
of this class of materials is concentrated on the applied-oriented
investigations, missing the systematic studies of their electronic structure
using different electron spectroscopy methods, which provide direct insight in
these properties. Such experiments can give a direct comparison with the
available theoretical data and link them to the available application-oriented
studies, providing the route to improve their properties 33.

Here, we present systematic x-ray photoelectron spectroscopy studies of the
representative TMT layered material, NiPS3. Using density functional theory
(DFT) calculations, it is shown that NiPS3 is an AFM semiconductor in its
ground state. It has a band gap of more than $2$ eV and the magnetic
properties of this material are described with a model Hamiltonian of the XY-
type. Near-edge x-ray absorption spectroscopy (NEXAFS) and resonant
photoelectron spectroscopy (ResPES) performed around the respective absorption
edges allow to discriminate between different atomic contributions in the
density of states (DOS) for valence and conduction bands. Our analysis of the
NEXAFS and ResPES data for the Ni $L_{2,3}$ absorption edge indicate that
NiPS3 can be described as a charge-transfer insulator according to the Zaanen-
Sawatzky-Allen scheme with strong correlation effects in the valence band of
this material. Details describing experimental and theoretical approaches used
in the present study as well as additional data are presented in the
Supplementary Information.

Layered NiPS3 bulk crystals in our study were synthesized using chemical vapor
transport (CVT) method from the stoichiometric amounts of elements and
obtained several-$\mathrm{mm}^{2}$ samples were characterized using different
bulk- and surface-sensitive techniques. Fig. 1(b,c) show x-ray diffraction
(XRD) plot, optical images and transmission electron microscopy (TEM) data for
NiPS3 crystal and these data are in perfect agreement with referenced and
previously published data 27, 34, 35 demonstrating high-quality of samples
studied in the present work (see also Fig. S1 and Fig. S2 in Supplementary
Information for the high-resolution cross-section TEM images and energy
dispersive x-ray (EDX) analysis data collected in the TEM measurements,
respectively). The distance between single layers of NiPS3 extracted from the
high-resolution TEM data is $6.357$ Å, which is in very good agreement with
the value of $6.301$ Å obtained in our DFT calculations for fully relaxed 3D
structure. [The complete crystallographic data for the fully optimized 3D
structure of NiPS3 are presented in Table S3.] Raman spectra of NiPS3 are
shown in Fig. S4 of Supplementary Information and also in very good agreement
with previously presented results, which were used for the assignment of the
observed scattering peaks 36, 37, 34, 38.

NiPS3 was characterized using different surface-sensitive methods and Fig.
1(d) presents low-energy electron diffraction (LEED) image of the freshly-
cleaved and UHV-degassed sample. It shows clear hexagonal diffraction spots,
which are due to the hexagonal unit cell of NiPS3 in the $xy$-plane
perpendicular to the impinge electron beam. The presented LEED image was
collected in the low-current electron-beam regime at the primary electron
energy of $95$ eV indicating high surface crystallographic quality with long-
range order. Our reference measurements using Au-crystal allow to conclude
that presented image shows 1st and 2nd diffraction orders for the NiPS3
$xy$-surface.

The results of x-ray photoelectron spectroscopy (XPS) characterization and
NEXAFS measurements at the Ni, S, and P $L_{2,3}$ absorption edges for NiPS3
are summarized in Fig. 2. Survey XPS spectrum measured for NiPS3 shows a
series of characteristic emission peaks (without any additional emission lines
from possible contaminants), which can be easily identified, and they are
marked in panel (a). The characteristic multiplet structure for the Ni $2p$
spectrum (Fig. 2(b)) is in very good agreement with the recent data for NiPS3
39 and NiGa2S4 40, where multiplet modeling using the NiS6 cluster gives the
dominant Ni $d^{9}\underline{L}$ ground state in both cases ($\underline{L}$
marks the hole on the ligand site). S $2p$ and P $2p$ XPS spectra (Fig.
2(c,d)) demonstrate clear single spin-orbit-slit doublets indicating strongly
localized covalent bonds between P and S atoms forming (P2S6)4- anions.

The Ni $L_{2,3}$ NEXAFS spectrum (Fig. 2(e) and Fig. 3(a)) demonstrates
characteristic absorption edge for two spin-orbit split components separated
by $\approx 17$ eV. According to the Zaanen-Sawatzky-Allen scheme for
transition-metal compounds (insulators and semiconductors) 41, their
electronic structure can be described with two parameters, the $d-d$
correlation energy ($U_{dd}$) and charge transfer energy between $d$-states of
metal and $p$-states of ligand ($\Delta$). Simplifying all considerations, one
can say that if $U_{dd}<\Delta$, then compound is the Mott-Hubbard type
insulator, while it is the charge-transfer type insulator if $U_{dd}>\Delta$.
For NiPS3 the Ni $L_{2,3}$ NEXAFS spectra were simulated in Ref. 39 and good
agreement with experimental data was suggested for $U_{dd}=5$ eV and
$\Delta=-1$ eV, indicating that this compound is the so-called self-doped
negative charge-transfer insulator. In our simulations (see computational
details and Fig. S5 in Supplementary Information) using $U_{dd}=6$ eV, the
best fit is achieved for charge transfer energy between $\Delta=0.5$ eV and
$\Delta=1$ eV giving simple charge-transfer insulating state for NiPS3.
Besides the correctly shown two-peaks satellite structure above the $L_{3}$
absorption edge, with these parameters our simulations also correctly
reproduce the high-energy shoulder for these line as well as all fine
structures in the NEXAFS spectra. Our DFT calculations with $U_{dd}=6$ eV lead
to the correct crystallographic structure of NiPS3 with simultaneous
adequately reproduced value for $T_{N}=148$ K (vs. $T_{N}=155$ K obtained in
the experiment 42). In this case, the respective value for the occupation
number for $d$ orbitals of the Ni2+ ion is $8.2$ and the calculated value
$1.6\,\mu_{B}$ for magnetic moment of Ni ions is closer to the formal value of
$2\,\mu_{B}$ compared to the previously reported results 39. Our DFT
calculations also give significantly wider band gap for NiPS3 – 2.19 eV (vs.
$\approx 1$ eV in Ref. 39) together with correctly calculated optical
absorption spectra (see Fig. S6 in Supplementary Information). Also our DFT
calculations give the different distribution of the PDOS in the valence band
compared to the previous results (see Fig. 1(a) and Fig. S7 for present DFT
results obtained using PBE$+U+$D2 and HSE06+D2 functionals, respectively). One
has to note that value $\Delta=0.95$ eV was recently obtained during correct
simulation of the resonant inelastic x-ray scattering spectra of NiPS3 17
supporting our results.

The S $L_{2,3}$ and P $L_{2,3}$ NEXAFS spectra of NiPS3 shown in Fig. 2(f) and
(g), respectively, have a very rich structure consisting of sharp peaks. These
spectra reflect the transitions from the ligand spin-orbit split $2p_{3/2}$
and $2p_{1/2}$ core levels onto the unoccupied $d$ and $s$ states. The
interpretation of these spectra is quite tedious and ambiguous procedure,
however one can conclude that the sharp double-peak structure in the S
$L_{2,3}$ ($163.6$ eV and $164.75$ eV) and P $L_{2,3}$ ($132.95$ eV and
$133.7$ eV) spectra reflects the electron transitions into the first
unoccupied hybrid $3s$-like antibonding state formed by S and P 43, 44. The
first peak in the S $L_{2,3}$ spectrum at $162.45$ eV and the respective
shoulder in the P $L_{2,3}$ spectrum at $132.2$ eV most probably correspond to
the states in the conduction band peaked at $E-E_{F}\approx 2.4$ eV. Further
peaks in the NEXAFS spectra observed at higher photon energies can be assigned
either to the so-called “echo” or shadow effect of the spin-orbit split
features due to multiple scattering or to electron transitions to a mixed-
valence band states 43, 44, 45.

The results of the ResPES measurements at the Ni $L_{2,3}$ absorption edge for
NiPS3 are presented in Fig. 3, supporting the description of this material as
the charge-transfer insulator: In panel (a) the reference Ni $L_{2,3}$ NEXAFS
spectrum is shown and in panel (b) a series of photoemission spectra taken at
the particular photon energies marked by the vertical arrows in (a) is
presented. Previous considerations of the transition metal oxides and sulfides
46, 47, 48, 49, 50 show that in case of the Mott-Hubbard insulator with
$U_{dd}<\Delta$ (e. g., Mn and Fe oxides) the core-level and valence band XPS
spectra consist of the low-binding-energy intense peak originated mainly from
the $3d^{n-1}$ final state and the satellite $3d^{n}\underline{L}$ structure
at larger binding energies. For the charge-transfer insulator with
$U_{dd}>\Delta$ (e. g., Ni oxide), the low-binding-energy intense peak is
formed by the configuration mixing of the $3d^{n}\underline{L}$ and
$3d^{n}\underline{Z}$ final states ($\underline{Z}$ is the bound states; e.
g., for NiO its a Zhang-Rice doublet bound states) 50 and the peak associated
with the $3d^{n-1}$ final state is located at larger binding energies with
smaller intensities. (Here one has to note that strong hybridization between
two final states does not allow to clearly separate them in the photoemission
spectra of the valence band region, requiring detailed fitting of the
experimental data using cluster model calculations.)

In the ResPES measurements at the Ni $L_{2,3}$ absorption edge for NiPS3 the
photoemission intensity is the result of interference of two emission channels
51: (i) a direct photoemission $2p^{6}3d^{n}+h\nu\rightarrow 2p^{6}3d^{n-1}+e$
and (ii) a photoabsorption process followed by a participator Coster-Kronig
decay $2p^{6}3d^{n}+h\nu\rightarrow 2p^{5}3d^{n+1}\rightarrow
2p^{6}3d^{n-1}+e$ and the initial and final states for these two photoemission
channels are identical. Also the interference between these excitation
channels leads to the Fano-type resonance for states with the $3d^{n-1}$
final-state character 51. For the $2p\rightarrow 3d$ ResPES experiments the
second emission channel (Coster-Kronig decay) prevails in the spectra,
therefore giving the position of the $3d^{n-1}$ peak relative to the position
of the $3d^{n}\underline{L}$ line.

The off-resonance spectrum of NiPS3 measured at $h\nu=843.34$ eV (spectrum $1$
in Fig. 3(b) and its intensity multiplied by factor 20 is shown by the dotted
line) demonstrates a broad two-peaks structure in the range of
$E-E_{VBM}=-5...0$ eV and additional satellite structure at larger binding
energies. If we assume that NiPS3 is a charge-transfer insulator with
$U_{dd}>\Delta$ (according to the previous considerations and considering the
similar ion state for charge-transfer insulating state of NiO), then the first
structure can be assigned to the $3d^{n}\underline{L}$ final state and other
peaks to the $3d^{n-1}$ final state. The on-resonance photoemission spectrum
measured with the photon energy $h\nu=851.94$ eV corresponding to the $L_{3}$
absorption edge (spectrum $3$ in Fig. 3(b)) demonstrates the drastic increase
of the photoemission intensity by factor of $20$ for the two-peaks structure
in the range $E-E_{VBM}=-5...0$ eV and by factor of $80$ for the peaks at
large binding energies. However, as it was discussed earlier, in the ResPES
experiments the $3d^{n-1}$ Coster-Kronig decay channel mainly contributes in
the photoemission spectra. Thus we can conclude that the huge peak at
$E-E_{VBM}\approx-6$ eV and the additional satellites structure at
$E-E_{VBM}\approx-12$ eV correspond to the $3d^{n-1}$ final state and confirm
that NiPS3 has to be described as the charge-transfer insulator material. The
similar consideration is valid for the on-resonance photoemission spectrum
measured at the $L_{2}$ absorption edge (spectrum $10$ in Fig. 3(b)). The
obtained results which allow to describe NiPS3 as charge-transfer insulator
are supported by the DFT results. In the DOS plot (Fig. 1(a)), the top of the
valence band is dominated by the S $p$ states with the Ni $3d$ states
spreading between $-3.5$ eV and $-6.5$ eV. At the same time the bottom of the
conduction band is formed mainly by unoccupied Ni $3d$ states, leading to the
description of NiPS3 as a charge-transfer insulator according to the Zaanen-
Sawatzky-Allen scheme (see Fig. S8) and confirming experimental spectroscopic
results. [Here we have to note that the correct simulation of the ResPES
spectra on the basis of the available DOS data is only possible in the many-
electron approximation, which might be a topic for further studies.] The
presented results are in a rather good agreement with the experimental and
calculated ResPES spectra of NiO (Ni2+) in the charge-transfer insulating
state, where similar spectral behavior was found 49, 50.

In summary, we performed systematic studies of the crystallographic structure
and electronic properties of layered van der Waals material NiPS3 using
different ex situ and in situ UHV surface science techniques, including LEED
and electron spectroscopy methods, like XPS, NEXAFS and ResPES. Our structural
investigations demonstrate very high quality of the studied crystals allowing
to obtain clear diffraction spots for the freshly-cleaved surface of
insulating NiPS3. Further results obtained using electron spectroscopy methods
indicate strong correlations effects in this material and we show that all
electron spectroscopy data can be successfully described using cluster based
approach for the Ni2+ ion in octahedral field. NEXAFS Ni $L_{2,3}$ data show
that NiPS3 can be described as a transfer-charge insulator with
$U_{dd}>\Delta$. This is also supported by the systematic $2p\rightarrow 3d$
ResPES data collected around the Ni $L_{2,3}$ absorption edge, which clearly
demonstrate the strong resonant behavior for the $3d^{n-1}$ photoemission
final state. As deduced from these data, this final state is located at large
binding energies compared to those for $3d^{n}\underline{L}$ final state
located at the top of the valence band. The present results are first
systematic electron spectroscopy data allowing to understand the electronic
structure of these new layered materials and are of paramount importance for
the further studies giving a strong basis for the understanding their
transport, optical, and catalytic properties.

This work was supported by the National Natural Science Foundation of China
(Grant No. 21973059). A.M. acknowledges the BMBF (grant No. 05K19KER). D.S.
acknowledges the BMBF (grant no. 0519ODR). A.T. acknowledges the Ministry of
Science and Higher Education of the Russian Federation no. 0852-2020-0019
(State assignment in the field of scientific activity, Southern Federal
University, 2020 project no. BAZ0110/20-1-03EH) and resources of the center of
the collective use of SFedU “High-Resolution Electron Microscopy”. We thank
HZB for the allocation of synchrotron radiation beamtime and for support
within the bilateral Russian-German Laboratory program. The North-German
Supercomputing Alliance (HLRN) is acknowledged for providing computer time.

The following files are available free of charge

  * •

Description of experimental methods, sample preparation and DFT calculations.
Additional experimental and theoretical data (PDF) can be downloaded via link:
https://pubs.acs.org/doi/10.1021/acs.jpclett.1c00394

## References

  * Geim 2009 Geim, A. Graphene: Status and Prospects. _Science_ 2009, _324_ , 1530–1534
  * Das Sarma et al. 2011 Das Sarma, S.; Adam, S.; Hwang, E. H.; Rossi, E. Electronic Transport in Two-Dimensional Graphene. _Rev. Mod. Phys._ 2011, _83_ , 407–470
  * Geim and Grigorieva 2014 Geim, A. K.; Grigorieva, I. V. Van der Waals Heterostructures . _Nature_ 2014, _499_ , 419–425
  * Zhuang et al. 2015 Zhuang, J.; Xu, X.; Feng, H.; Li, Z.; Wang, X.; Du, Y. Honeycomb Silicon: A Review of Silicene. _Sci. Bull._ 2015, _60_ , 1551–1562
  * Carvalho et al. 2016 Carvalho, A.; Wang, M.; Zhu, X.; Rodin, A. S.; Su, H.; Castro Neto, A. H. Phosphorene: From Theory to Applications. _Nat. Rev. Mater._ 2016, _1_ , 16061
  * Manzeli et al. 2017 Manzeli, S.; Ovchinnikov, D.; Pasquier, D.; Yazyev, O. V.; Kis, A. 2D Transition Metal Dichalcogenides. _Nat. Rev. Mater._ 2017, _2_ , 17033
  * Bae et al. 2010 Bae, S.; Kim, H.; Lee, Y.; Xu, X.; Park, J.-S.; Zheng, Y.; Balakrishnan, J.; Lei, T.; Kim, H. R.; Song, Y. I. et al. Roll-to-Roll Production of 30-inch Graphene Films for Transparent Electrodes. _Nat. Nanotech._ 2010, _5_ , 574–578
  * Chhowalla et al. 2013 Chhowalla, M.; Shin, H. S.; Eda, G.; Li, L.-J.; Loh, K. P.; Zhang, H. The Chemistry of Two-Dimensional Layered Transition Metal Dichalcogenide Nanosheets. _Nat. Chem._ 2013, _5_ , 263–275
  * Ryu et al. 2014 Ryu, J.; Kim, Y.; Won, D.; Kim, N.; Park, J. S.; Lee, E.-K.; Cho, D.; Cho, S.-P.; Kim, S. J.; Ryu, G. H. et al. Fast Synthesis of High-Performance Graphene Films by Hydrogen-Free Rapid Thermal Chemical Vapor Deposition. _ACS Nano_ 2014, _8_ , 950–956
  * Dedkov and Voloshina 2015 Dedkov, Y.; Voloshina, E. Graphene Growth and Properties on Metal Substrates. _J. Phys.: Condens. Matter_ 2015, _27_ , 303002
  * Li et al. 2018 Li, H.; Li, Y.; Aljarb, A.; Shi, Y.; Li, L.-J. Epitaxial Growth of Two-Dimensional Layered Transition-Metal Dichalcogenides: Growth Mechanism, Controllability, and Scalability. _Chem. Rev._ 2018, _118_ , 6134–6150
  * Dedkov and Voloshina 2020 Dedkov, Y.; Voloshina, E. Epitaxial Graphene/Ge Interfaces: A Minireview. _Nanoscale_ 2020, _12_ , 11416–11426
  * Cao et al. 2018 Cao, Y.; Fatemi, V.; Fang, S.; Watanabe, K.; Taniguchi, T.; Kaxiras, E.; Jarillo-Herrero, P. Unconventional Superconductivity in Magic-Angle Graphene Superlattices. _Nat. Rev. Mater._ 2018, _556_ , 43
  * Andrei and MacDonald 2020 Andrei, E. Y.; MacDonald, A. H. Graphene Bilayers with a Twist. _Nat. Mater._ 2020, _19_ , 1265–1275
  * Xi 2015 Xi, X. Strongly Enhanced Charge-Density-Wave Order in Monolayer NbSe2. _Nat. Nanotechnol._ 2015, _10_ , 765
  * Hossain et al. 2017 Hossain, M.; Zhao, Z.; Wen, W.; Wang, X.; Wu, J.; Xie, L. Recent Advances in Two-Dimensional Materials with Charge Density Waves: Synthesis, Characterization and Applications. _Crystals_ 2017, _7_ , 298–19
  * Kang et al. 2020 Kang, S.; Kim, K.; Kim, B. H.; Kim, J.; Sim, K. I.; Lee, J.-U.; Lee, S.; Park, K.; Yun, S.; Kim, T. et al. Coherent Many-Body Exciton in Van der Waals Antiferromagnet NiPS3. _Nature_ 2020, _583_ , 785–789
  * Lee et al. 2016 Lee, J.-U.; Lee, S.; Ryoo, J. H.; Kang, S.; Kim, T. Y.; Kim, P.; Park, C.-H.; Park, J.-G.; Cheong, H. Ising-Type Magnetic Ordering in Atomically Thin FePS3. _Nano Lett._ 2016, _16_ , 7433–7438
  * Wang et al. 2016 Wang, X.; Du, K.; Liu, Y. Y. F.; Hu, P.; Zhang, J.; Zhang, Q.; Owen, M. H. S.; Lu, X.; Gan, C. K.; Sengupta, P. et al. Raman Spectroscopy of Atomically Thin Two-Dimensional Magnetic Iron Phosphorus Trisulfide (FePS3) Crystals. _2D Mater._ 2016, _3_ , 031009
  * Huang et al. 2017 Huang, B.; Clark, G.; Navarro-Moratalla, E.; Klein, D. R.; Cheng, R.; Seyler, K. L.; Zhong, D.; Schmidgall, E.; McGuire, M. A.; Cobden, D. H. et al. Layer-Dependent Ferromagnetism in a Van der Waals Crystal Down to the Monolayer Limit. _Nature_ 2017, _546_ , 270
  * Gong et al. 2017 Gong, C.; Li, L.; Li, Z.; Ji, H.; Stern, A.; Xia, Y.; Cao, T.; Bao, W.; Wang, C.; Wang, Y. et al. Discovery of Intrinsic Ferromagnetism in Two-Dimensional Van der Waals Crystals. _Nature_ 2017, _546_ , 265
  * Deng et al. 2018 Deng, Y.; Yu, Y.; Song, Y.; Zhang, J.; Wang, N. Z.; Sun, Z.; Yi, Y.; Wu, Y. Z.; Wu, S.; Zhu, J. et al. Gate-Tunable Room-Temperature Ferromagnetism in Two-Dimensional Fe3GeTe2. _Nature_ 2018, _563_ , 94
  * O’Hara et al. 2018 O’Hara, D. J.; Zhu, T.; Trout, A. H.; Ahmed, A. S.; Luo, Y. K.; Lee, C. H.; Brenner, M. R.; Rajan, S.; Gupta, J. A.; McComb, D. W. et al. Room Temperature Intrinsic Ferromagnetism in Epitaxial Manganese Selenide Films in the Monolayer Limit. _Nano Lett._ 2018, _18_ , 3125–3131
  * Susner et al. 2017 Susner, M. A.; Chyasnavichyus, M.; McGuire, M. A.; Ganesh, P.; Maksymovych, P. Metal Thio- and Selenophosphates as Multifunctional Van der Waals Layered Materials. _Adv. Mater._ 2017, _29_ , 1602852
  * Wang et al. 2018 Wang, F.; Shifa, T. A.; Yu, P.; He, P.; Liu, Y.; Wang, F.; Wang, Z.; Zhan, X.; Lou, X.; Xia, F. et al. New Frontiers on van der Waals Layered Metal Phosphorous Trichalcogenides. _Adv. Funct. Mater._ 2018, _28_ , 1802151
  * Zhu et al. 2020 Zhu, M.; Kou, H.; Wang, K.; Wu, H.; Ding, D.; Zhou, G.; Ding, S. Promising Functional Two-Dimensional Lamellar Metal Thiophosphates: Synthesis Strategies, Properties and Applications. _Mater. Horiz._ 2020, _529_ , 185
  * Du et al. 2016 Du, K.-z.; Wang, X.-z.; Liu, Y.; Hu, P.; Utama, M. I. B.; Gan, C. K.; Xiong, Q.; Kloc, C. Weak Van der Waals Stacking, Wide-Range Band Gap, and Raman Study on Ultrathin Layers of Metal Phosphorus Trichalcogenides. _ACS Nano_ 2016, _10_ , 1738–1743
  * Yang et al. 2020 Yang, J.; Zhou, Y.; Guo, Q.; Dedkov, Y.; Voloshina, E. Electronic, Magnetic and Optical Properties of MnPX3 (X = S, Se) Monolayers with and without Chalcogen Defects: A First-Principles Study. _RSC Adv._ 2020, _10_ , 851–864
  * Sivadas et al. 2015 Sivadas, N.; Daniels, M. W.; Swendsen, R. H.; Okamoto, S.; Xiao, D. Magnetic Ground State of Semiconducting Transition-Metal Trichalcogenide Monolayers. _Phys. Rev. B_ 2015, _91_ , 235425
  * Chittari et al. 2016 Chittari, B. L.; Park, Y.; Lee, D.; Han, M.; MacDonald, A. H.; Hwang, E.; Jung, J. Electronic and Magnetic Properties of Single-Layer MPX3 Metal Phosphorous Trichalcogenides. _Phys. Rev. B_ 2016, _94_ , 184428
  * Li et al. 2014 Li, X.; Wu, X.; Yang, J. Half-Metallicity in MnPSe3 Exfoliated Nanosheet with Carrier Doping. _J. Am. Chem. Soc._ 2014, _136_ , 11065–11069
  * Zhang et al. 2016 Zhang, X.; Zhao, X.; Wu, D.; Jing, Y.; Zhou, Z. MnPSe3 Monolayer: A Promising 2D Visible-Light Photohydrolytic Catalyst with High Carrier Mobility. _Adv. Sci._ 2016, _3_ , 1600062
  * Yang et al. 2020 Yang, J.; Zhou, Y.; Dedkov, Y.; Voloshina, E. Dirac Fermions in Half-Metallic Ferromagnetic Mixed Cr1-xMxPSe3 Monolayers. _Adv. Theory Simul._ 2020, _3_ , 2000228
  * Jenjeti et al. 2018 Jenjeti, R. N.; Kumar, R.; Austeria, M. P.; Sampath, S. Field Effect Transistor Based on Layered NiPS3. _Sci. Rep._ 2018, _8_ , 8586
  * Xue et al. 2018 Xue, S.; Chen, L.; Liu, Z.; Cheng, H.-M.; Ren, W. NiPS3 Nanosheet–Graphene Composites as Highly Efficient Electrocatalysts for Oxygen Evolution Reaction. _ACS Nano_ 2018, _12_ , 5297–5305
  * Bernasconi et al. 1988 Bernasconi, M.; Marra, G. L.; Benedek, G.; Miglio, L.; Jouanne, M.; Julien, C.; Scagliotti, M.; Balkanski, M. Lattice Dynamics of Layered MPX3 (M = Mn, Fe, Ni, Zn; X = S, Se) Compounds. _Phys. Rev. B_ 1988, _38_ , 12089–12099
  * Kuo et al. 2016 Kuo, C.-T.; Neumann, M.; Balamurugan, K.; Park, H. J.; Kang, S.; Shiu, H. W.; Kang, J. H.; Hong, B. H.; Han, M.; Noh, T. W. et al. Exfoliation and Raman Spectroscopic Fingerprint of Few-Layer NiPS3 Van der Waals Crystals. _Sci. Rep._ 2016, _6_ , 20904
  * Kim et al. 2019 Kim, K.; Lim, S. Y.; Lee, J.-U.; Lee, S.; Kim, T. Y.; Park, K.; Jeon, G. S.; Park, C.-H.; Park, J.-G.; Cheong, H. Suppression of Magnetic Ordering in XXZ-Type Antiferromagnetic Monolayer NiPS3. _Nat. Commun._ 2019, _10_ , 345
  * Kim et al. 2018 Kim, S. Y.; Kim, T. Y.; Sandilands, L. J.; Sinn, S.; Lee, M.-C.; Son, J.; Lee, S.; Choi, K.-Y.; Kim, W.; Park, B.-G. et al. Charge-Spin Correlation in Van der Waals Antiferromagnet NiPS3. _Phys. Rev. Lett._ 2018, _120_ , 136402
  * Takubo et al. 2007 Takubo, K.; Mizokawa, T.; Son, J. Y.; Nambu, Y.; Nakatsuji, S.; Maeno, Y. Unusual Superexchange Pathways in an NiS2 Triangular Lattice with Negative Charge-Transfer Energy. _Phys. Rev. Lett._ 2007, _99_ , 037203
  * Zaanen et al. 1985 Zaanen, J.; Sawatzky, G. A.; Allen, J. W. Band Gaps and Electronic Structure of Transition-Metal Compounds. _Phys. Rev. Lett._ 1985, _55_ , 418–421
  * Le Flem et al. 1982 Le Flem, G.; Brec, R.; Ouvard, G.; Louisy, A.; Segransan, P. Magnetic Interactions in the Layer Compounds MPX3 (M = Mn, Fe, Ni; X = S, Se). _J. Phys. Chem. Solids_ 1982, _43_ , 455–461
  * Farrell et al. 2002 Farrell, S. P.; Fleet, M. E.; Stekhin, I. E.; Kravtsova, A.; Soldatov, A. V.; Liu, X. Evolution of Local Electronic Structure in Alabandite and Niningerite Solid Solutions [(Mn,Fe)S, (Mg,Mn)S, (Mg,Fe)S] Using Sulfur K- and L-edge XANES Spectroscopy. _Am. Mineral._ 2002, _87_ , 1321–1332
  * Kruse et al. 2009 Kruse, J.; Leinweber, P.; Eckhardt, K.-U.; Godlinski, F.; Hu, Y.; Zuin, L. Phosphorus L2,3-Edge XANES: Overview of Reference Compounds. _J. Synchrotron Rad._ 2009, _16_ , 247–259
  * Yang et al. 2012 Yang, S.; Wang, D.; Liang, G.; Yiu, Y. M.; Wang, J.; Liu, L.; Sun, X.; Sham, T.-K. Soft X-ray XANES Studies of Various Phases Related to LiFePO4 Based Cathode Materials. _Energy Environ. Sci._ 2012, _5_ , 7007
  * Fujimori and Minami 1984 Fujimori, A.; Minami, F. Valence-Band Photoemission and Optical Absorption in Nickel Compounds. _Phys. Rev. B_ 1984, _30_ , 957–971
  * Fujimori et al. 1990 Fujimori, A.; Kimizuka, N.; Akahane, T.; Chiba, T.; Kimura, S.; Minami, F.; Siratori, K.; Taniguchi, M.; Ogawa, S.; Suga, S. Electronic Structure of MnO. _Phys. Rev. B_ 1990, _42_ , 7580–7586
  * Okada et al. 1992 Okada, K.; Kotani, A.; Thole, B. T. Charge Transfer Satellites and Multiplet Splitting in X-ray Photoemission Spectra of Late Transition Metal Halides. _J. El. Spectr. Rel. Phenom._ 1992, _58_ , 325–343
  * Tjernberg et al. 1996 Tjernberg, O.; Söderholm, S.; Karlsson, U. O.; Chiaia, G.; Qvarford, M.; Nylén, H.; Lindau, I. Resonant Photoelectron Spectroscopy on NiO. _Phys. Rev. B_ 1996, _53_ , 10372–10376
  * Taguchi et al. 2008 Taguchi, M.; Matsunami, M.; Ishida, Y.; Eguchi, R.; Chainani, A.; Takata, Y.; Yabashi, M.; Tamasaku, K.; Nishino, Y.; Ishikawa, T. et al. Revisiting the Valence-Band and Core-Level Photoemission Spectra of NiO. _Phys. Rev. Lett._ 2008, _100_ , 206401
  * Davis and Feldkamp 1981 Davis, L. C.; Feldkamp, L. A. Resonant Photoemission Involving Super-Coster-Kronig Transitions. _Phys. Rev. B_ 1981, _23_ , 6239–6253

Figure 1: (a) Top and side views of the crystallographic structure of the
NiPS3 single layer. Atom-projected partial DOS for bulk NiPS3 in the AFM
ground state is shown on the left-hand side of the panel. (b) XRD patterns of
NiPS3 measured at room temperature. Insets show photo of the sample and its
optical microscopy image, respectively. (c) TEM image of the NiPS3. Inset
shows the respective large scale TEM image of several-$\mu\mathrm{m}^{2}$
NiPS3 sheet. (d) MCP-LEED image of the NiPS3(001) surface acquired at $95$ eV.

Figure 2: XPS spectra of NiPS3 collected at $h\nu=1100$ eV: (a) survey, (b) Ni
$2p$, (c) S $2p$, (d) P $2p$. NEXAFS spectra of NiPS3 collected in the TEY
mode: (e) Ni $L_{2,3}$, (f) S $L_{2,3}$, (g) P $L_{2,3}$.

Figure 3: Results of ResPES experiments for NiPS3: (a) the reference Ni
$L_{2,3}$ NEXAFS spectrum and (b) a series of photoemission spectra taken at
the particular photon energies marked by the corresponding vertical arrows in
panel (a). Spectrum measured at the off-resonance photon energy (spectrum 1)
with intensity multiplied by factor 20 is also shown by the doted line. Solid
thin curves drawn on top of the spectra $3-12$ and marked by a1 and a2 note
the respective Auger-decay lines. All ResPES spectra are shifted in vertical
direction for clarity. Zero of the energy scale for ResPES data is referenced
to the valence band maximum edge for NiPS3 extracted from the off-resonance
spectra.


# Non-parametric Bayesian mixture model to study adverse events of COVID-19
vaccines

Ali Turfah
Department of Biostatistics
University of Michigan
Ann Arbor, MI 48109
<EMAIL_ADDRESS>
&Xiaoquan Wen
Department of Biostatistics
University of Michigan
Ann Arbor, MI 48109
<EMAIL_ADDRESS>
&Lili Zhao∗
Research Institute
Corewell Health’s Beaumont Hospital
Royal Oak, MI 48073
<EMAIL_ADDRESS>
Both Xiaoquan Wen and Lili Zhao are corresponding authors.

(March 24, 2023)

###### Abstract

The vaccine adverse event reporting system (VAERS) is a vital resource for
post-licensure vaccine safety monitoring and has played a key role in
assessing the safety of COVID-19 vaccines. However it is difficult to properly
identify rare adverse events (AEs) associated with vaccines due to small or
zero counts. We propose a Bayesian model with a Dirichlet Process Mixture
prior to improve accuracy of the AE estimates with small counts by allowing
data-guided information sharing between AE estimates. We also propose a
negative control procedure embedded in our Bayesian model to mitigate the
reporting bias due to the heightened awareness of COVID-19 vaccines, and use
it to identify associated AEs as well as associated AE groups defined by the
organ system in the Medical Dictionary for Regulatory Activities (MedDRA)
ontology. The proposed model is evaluated using simulation studies, in which
it outperforms baseline models without information sharing and is applied to
study the safety of COVID-19 vaccines using VAERS data.

_K_ eywords Vaccine safety $\cdot$ Dirichlet process mixture $\cdot$ Adverse
event $\cdot$ COVID-19 vaccines $\cdot$ VAERS $\cdot$ Enrichment test

##  1 Background

As of January 1, 2023, COVID-19 vaccines have been administered to over 600
million people in the United States [1]. However, concern about the safety of
the vaccines is among the primary factors contributing to vaccine hesitancy
[2]. There are three COVID-19 vaccines approved by the US Food and Drug
Administration (FDA) under emergency use authorization (EUA) in the United
States: two mRNA-based vaccines—the BNT162b2 (Pfizer–BioNTech) and the
mRNA-1273 (Moderna) vaccines—and the Ad26.COV2.S (Johnson & Johnson–Janssen)
vaccine.

Due to limited sample size, restrictive inclusion criteria, and the limited
duration of follow-up in phase-III trials, a chief concern is that rare or
serious adverse events associated with the COVID-19 vaccine have not yet been
identified. In addition, despite being theoretically safe, the two mRNA-based
vaccines (BNT162b2 and mRNA-1273) are coming to the market for the first time;
their safety can only be verified by long-term longitudinal surveillance.

The Centers for Disease Control and Prevention (CDC) and the U.S. Food and
Drug Administration (FDA) conduct post-licensure vaccine safety monitoring
using the Vaccine Adverse Event Reporting System (VAERS) [3]. VAERS accepts
spontaneous reports of suspected vaccine adverse events after administration
of any vaccine licensed in the United States from 1990 to the present. As a
national public health surveillance resource, VAERS is key for ensuring the
safety of vaccines. We define an adverse event (AE) signal as a potential
vaccine safety problem, indicating that an AE might be caused by vaccination
and warrants further investigation or action. Several AE signals for COVID-19
vaccines have been identified using VAERS data, including anaphylaxis and
myocarditis after receipt of mRNA-based COVID-19 vaccines [4, 5], and
Guillain-Barré syndrome and Cerebral Venous Sinus Thrombosis (CVST) after
receipt of the Ad26.COV2.S vaccine [4, 6].

As a passive surveillance system, VAERS lacks information on the total number
of subjects vaccinated, the total number of vaccine recipients who experience
an AE, or the incidence of AEs in unvaccinated subjects. However, the
proportion of reports mentioning a specific AE for the COVID-19 vaccines can
be compared to the proportion of reports mentioning the same AE for a
reference vaccine—or all other vaccines (i.e., non-COVID-19 vaccines). This
type of approach is called disproportionality analysis. VAERS data is commonly
represented by a large contingency table for the counts of all vaccine and AE
pairs over a specified time period. Each row of this table represents a
vaccine and each column represents an AE. Each cell in the table contains the
number of VAERS reports that mention both the vaccine and the AE during the
time period. This table can be very large, with hundreds of rows for vaccines
and thousands of columns for AEs. To study the relationship between COVID-19
vaccines and each AE, the large table can be collapsed into a $2\times 2$
table with rows corresponding to reciept of a COVID-19 vaccine (Yes/No) and
two columns for the occurrence of the AE (Yes/No). The strength of the
vaccine-AE association can be measured using a ratio, such as a proportional
reporting ratio (PRR) or an reporting odds ratio (ROR) [7]. A larger ratio
($>1$) corresponds to a stronger association, which might indicate a vaccine
safety problem.

Many disproportionality methods have been used to conduct safety studies based
on the VAERS database [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22]. Generally these methods investigate each AE separately, assuming that the
AE counts follow a Poisson or Binomial distribution. As we know, AEs from
vaccines are often rare with sparse count data. In these separate analyses,
the ratio parameter for the rare AEs either cannot be estimated or will take
an extreme value with a large degree of uncertainty. Therefore, efficient
approaches are needed to improve estimation for rare AEs. The Gamma-Poisson
Shrinker (GPS) model [8] enhances the simple use of the separate AE analyses
by allowing information borrowing between similar rate estimates, thereby
reducing the sampling variation in the data and improving accuracy for the
rare AE estimation. GPS assumes that report counts for the AEs follows a
Poisson distribution with a rate representing the strength of the association
(i.e., rate $>1$ for a signaled AE). This approach uses a mixture of two
distributions to model the rates of “signaled" and “non-signaled" AEs,
assuming that majority of the rates are clustered at or below one (“non-
signaled" AEs) while the remaining rates are clustered at a rate above one
(“signaled" AEs). Although theoretically sound, it is practically restrictive
to approximate the distribution of the rates using only two distributions. The
Multi-item Gamma-Poisson Shrinker (MGPS) [16] extends the GPS model to adjust
for patient risk factors. However, MGPS can only handle categorical risk
factors and has increased difficulty estimating the rates for rare AEs because
it creates separate contingency tables for each risk factor level combination,
exacerbating the sparse counts already present in VAERS.

In this paper, we estimate the ratio parameter in a regression framework,
allowing for both categorical and continuous risk factors. We approximate the
distribution of the ratio parameters by a mixture of distributions where the
number of mixing components is automatically determined from the data.
Specifically, we assume that AE count data follow a binomial distribution,
modeling the AE count as a function of the vaccine group (e.g., COVID-19
vaccines Yes/No) and risk factors (such as age and gender) using Bayesian
logistic regression. We assign a Dirichlet process mixture (DPM) [23] prior to
the log odds ratio parameter; we refer to it as the log reporting odds ratio
(logROR) parameter in the VAERS data analysis. The DPM prior uses a mixture of
normal distributions to estimate logROR without determining the number of
mixing components beforehand; it efficiently accommodates signalled and non-
signalled AEs and allows for a more flexible distribution of ROR parameters
with the possibility of skewness and multimodality. The logROR parameters with
similar magnitudes (within a cluster) will be shrunk towards each other,
thereby borrowing information from each other to improve the accuracy of
logROR estimates, increasing the power to detect true signals while reducing
false positives.

Like all spontaneous reporting systems, VAERS data is vulnerable to reporting
bias [7]. Due to the heightened awareness around the COVID-19 vaccines, AEs
observed after administration of COVID-19 vaccines are more likely to be
reported than those from other vaccines, which will bias association estimates
and result in false positive associations. In order to mitigate this bias,
negative controls (NCs)—AEs that are known not to be causally related to any
vaccine [24, 25]—are often used, requiring the ROR estimates for AE signals to
be significantly larger than the ROR estimates of the NCs. Moreover, based on
MedDRA (Medical Dictionary for Regulatory Activities) [26], the largest AE
ontology resource describing disease relationships, we conduct enrichment
analysis to identify enriched AE groups where signaled AEs occur more
frequently compared to other groups. The proposed methods have been
implemented in JAGS and the source code is included in Supplementary Document
A.

The rest of the article is organized as follows. In Section 2 we develop the
DPM model. In Section 3 we present the results of our simulation studies. In
Section 4 we apply the proposed methods to study the safety of COVID-19
vaccines using the VAERS database. Concluding remarks are given in Section 5.

##  2 Methods

###  2.1 Bayesian Logistic Regression

Suppose there are a total of $N$ reports in VAERS and $J$ AEs mentioned in
these reports. Let $V_{i}\in\\{0,1\\}$ indicate the type of the vaccine (1 for
COVID-19 vaccines and 0 for non-COVID-19 vaccines) for report $i$
$(i=1,\ldots,N)$. Let $A_{ij}$ be the indicator of mentioning AE $j$
$(j=1,\ldots,J)$ in report $i$, i.e., $A_{ij}=1$ if AE $j$ is mentioned in
report $i$ and 0 otherwise. Let $\textbf{X}_{i}$ be a vector of covariates for
report $i$ (including an intercept term). We use a logistic regression to
estimate the strength of association between COVID-19 vaccines and AE $j$,

$logit(Pr(A_{ij}=1|\mathbf{X}_{i},V_{i}))=\boldsymbol{\alpha}_{j}^{\intercal}\mathbf{X}_{i}+\beta_{j}V_{i},~{}~{}~{}i=1,\cdots,N,$

where $\beta_{j}$ is the log reporting odds ratio (logROR) of AE $j$ when
comparing a COVID-19 vaccine to a non-COVID-19 vaccine, which is the parameter
of interest.

To allow information borrowing for estimating the logROR, we assign
$\beta_{j}$ a Dirichlet Process Mixture (DPM) of normals prior. That is,

$\displaystyle\beta_{1},\cdots,\beta_{J}$ $\displaystyle\sim G$ $\displaystyle
G$ $\displaystyle\sim DP(\gamma,G_{0}),$

where $G_{0}$ corresponds to a best guess for $G$ as a priori and $\gamma$
expresses confidence in this guess. The stick-breaking representation [27]
implies that $G\sim DP(\gamma,G_{0})$ is equivalent to

$\displaystyle G$
$\displaystyle=\sum_{k=1}^{\infty}\pi_{k}\textbf{N}(\mu_{k},\tau_{k})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(\mu_{k},\tau_{k})\sim G_{0}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\sum_{h=1}^{\infty}\pi_{h}=1$

The stick-breaking approach parametrizes the component probabilities $\pi_{k}$
as the proportion of remaining probability after considering
$\sum_{h=1}^{k-1}\pi_{h}$. The proportion of the remaining probabilities for
each component is assigned a $Beta(a_{0},b_{0})$ prior. Recent research has
focused on using the constructive definition of the Dirichlet process to
produce practical MCMC algorithms [28]. The principle is to approximate the
full process by truncating the DPM at a maximum number of components $K$, so
that

$G=\sum_{k=1}^{K}\pi_{k}\textbf{N}(\mu_{k},\tau_{k})$

While a large $K$ provides an accurate approximation to the full DPM, it
computationally expensive to fit. Strategies have been proposed to specify $K$
[23, 28]. We found that $K=5$ worked very well for the VAERS data in our
applications.

These $K$ components can be interpreted as clusters from which the $\beta_{j}$
of the AEs are drawn [23]. These clusters enables us to shrink very extreme
estimates of $\beta_{j}$ to more reasonable values as well as obtain more
stable estimates for AEs observed in a very small number of cases, neither of
which would be feasible if the AEs were to be modeled separately.

We then specify the hyperpriors for $\mu_{k}$ and $\tau_{k}$ in the $k^{th}$
component of the DPM based on [23, 29],

$\displaystyle\mu_{k}$ $\displaystyle\sim N(\mu_{base},\tau_{base})$
$\displaystyle\mu_{base}$ $\displaystyle\sim N(0,\tau_{0})$
$\displaystyle\tau_{base}$ $\displaystyle=\frac{1}{u^{2}}$ $\displaystyle u$
$\displaystyle\sim Uniform(0,f_{0})$ $\displaystyle\tau_{k}$
$\displaystyle\sim Gamma(r_{\tau},\lambda)$ $\displaystyle\lambda$
$\displaystyle\sim Gamma(r_{\lambda},\lambda_{0}),$

where $\tau_{0},f_{0},r_{\tau},r_{\lambda},\lambda_{0}$ are hyperparameters
set a priori. For other coefficients $\boldsymbol{\alpha}_{j}$ in the logistic
regression model, we use normal priors with mean 0 and precision
$\tau_{\alpha}$.

###  2.2 Bayesian Negative Control Approach to Identify AE Signals

In order to mitigate the reporting bias from analyses using spontaneous
reporting systems such as VAERS, we identify $36$ AE terms in VAERS as
negative controls (NCs). The list of NCs used in our analyses can be found int
Supplementary Document A. All of the NCs are known not to be causally related
to any vaccine and are reported at least 300 times in our dataset. The chosen
NCs occur between 300-16,300 times in our dataset. If there is no reporting
bias, the RORs of the NCs should be centered around one. With the
aforementioned over-reporting issue, RORs for these NCs will be centered
around a value larger than one. Therefore, we can not identify signal AEs
based solely on whether their ROR is larger than one or not. In [25] the
authors use ROR estimates from the NCs to model the null distribution and use
this empirically estimated distribution to yield properly calibrated p-values.

Our Bayesian model can easily incorporate these NCs to make proper posterior
inference. At each MCMC iteration, we obtain the logROR estimate for AE $j$
($\beta_{j}$) and compare it to the logROR estimates for all NC AEs in that
iteration ($\beta_{NC_{1}},\cdots,\beta_{NC_{36}}$). If $\beta_{j}$ is larger
than majority of the NC logROR estimates ($>35$ NCs is used in this paper)
then it is assigned a signal indicator of 1 for the iteration. Otherwise it is
assigned a signal indicator of 0 for the iteration. Across all iterations, We
compute the posterior probability that AE $j$ is a signal based on the
proportion of iterations in which the AE was assigned a signal indicator of 1.
Finally, we define an AE as a signal if this NC-adjusted posterior probability
is larger than a threshold—a 90% cutoff is used in this paper.

###  2.3 Bayesian Enrichment Test to Identify Enriched AE Groups

The MedDRA ontology [26] groups AEs into System Organ Classes (SOC) based on
aetiology, manifestation site, or purpose. For instance, “Dyspnea", “Atrial
Thrombosis", and “Myocardial Oedema" are all members of the “Cardiac Disorder"
SOC group. In addition to identifying individual AEs associated with the
COVID-19 vaccine, we are interested in identifying SOC groups which are
associated with the vaccine; for example, it is important to know that
COVID-19 vaccines are associated with cardiac disorders.

Let $\mathcal{S}$ denote the set of all SOC groups. To conduct the enrichment
analysis for a particular SOC group, $s$, a common approach is to create a
$2\times 2$ table similar to the one presented in Table 1. The elements of
this table are the counts of the $J$ AEs grouped based on their membership in
group $s$ as well as their signal AE status (determined based on a
significance test). If there are significantly more signaled AEs in $s$, based
on a Chi-square or Fisher’s exact test, then group $s$ is enriched.

SOC Group # AEs | Signal | Non-Signal | Total
---|---|---|---
$s$ (e.g. Cardiac Disorders) | $a$ | $J_{s}-a$ | $J_{s}$
$\mathcal{S}-s$ (All others) | $c$ | $J_{s^{*}}-c$ | $J_{s^{*}}$
Total | $a+c$ | $J-a-c$ | $J$
Table 1: Contingency table for SOC group enrichment analysis

In this paper, we propose an enrichment analysis in a Bayesian framework. We
construct this $2\times 2$ table in every MCMC iteration and perform the
enrichment analysis based on the resulting posterior distribution.
Specifically, at each MCMC iteration we obtain the signal indicators for AEs
within group $s$ ($a$ signals from $J_{s}$ AEs) as well as for AEs from all
other groups ($c$ signals from $J_{s^{*}}$ AEs). We then compute an odds ratio
(OR) based on Table 1. This OR represents an enrichment score (we refer to it
as the enrichment OR, or EOR); a larger EOR ($>1$) indicates stronger evidence
that $s$ is enriched. Across all MCMC iterations, we obtain the posterior
distribution of the EOR for $s$. Group $s$ is considered enriched if the lower
bound of the 95% posterior credible interval for the EOR is greater than 1 and
the posterior mean of the EOR is greater than a pre-specified value to ensure
clinical relevance.

Compared to the conventional approach described above, our method takes into
account the uncertainty in the measurement of signal significance of the AEs;
that is, signal indicators may vary across MCMC iterations while the
conventional approach has the signal indicators fixed in the enrichment
analysis. Moreover, it is very simple to implement in the Bayesian framework
since the test is constructed using by-products of posterior draws.

##  3 Simulation Study

| $\sigma=0.5$ |  | $\sigma=0.8$ |  | $\sigma=1$
---|---|---|---|---|---
| DIC | MSE | Coverage |  | DIC | MSE | Coverage |  | DIC | MSE | Coverage
DPM | 3,653 | 0.87 | 95.3% |  | 3,624 | 1.16 | 93.4% |  | 3,673 | 0.97 | 96.6%
IL-Informative | 3,669 | 0.93 | 95.7% |  | 3,649 | 1.21 | 94.7% |  | 3,675 | 1.05 | 95.0%
IL-Vague | 3,691 | 6.43 | 94.7% |  | 3,652 | 6.64 | 94.7% |  | 3,694 | 6.29 | 96.6%
Table 2: DPM and IL model performance in simulation studies.

We set up simulation studies to compare our proposed DPM model (allowing
information sharing between logROR parameters) against a logistic regression
with independent normal priors for the logROR parameters (referred to as the
IL model), which does not allow information sharing. We investigated two IL
models: one has a more informative prior on the logROR parameters (IL-
Informative), with precision of 0.1, whereas the other one (IL-Vague) has a
smaller precision of 0.01. To mimic the real data and have a distribution of
AE counts similar to that from the VAERS data, we took a random sample of 300
intercept terms from the mRNA (target) vs FLU (control) analysis (Section 4.2)
using a logistic regression model. We then simulated the logROR parameters
from a mixture of three normal distributions centered at -2.5, 0, and 2.5,
with size 50, 150 and 100, respectively. This corresponds to the belief that a
majority of the AEs will not be associated with any vaccine, however more AEs
will be associated with the target vaccine than with the control vaccine. We
also assumed the within-cluster standard deviation to be the same for the
three normal distributions, taking $\sigma=0.5,0.8,$ or $1$. As $\sigma$
increases, the clusters become less distinct (as can be seen from the heatmaps
in Figures 2(a), 2(b) and 2(c)) and so we expect the clustering ability of our
model to worsen.

We generated 250 datasets for each value of $\sigma$ under the above set-up.
All simulated datasets have sparse AE counts; approximately 33% of the AEs
have 0 counts for the control vaccine group. Across the three AE clusters
(centered at -2.5, 0, 2.5) between 32-33%, 11-12%, and approximately 4.5% of
the AEs have 0 counts for the target vaccine group respectively.

The precision for $\boldsymbol{\alpha}_{j}$ in both the DPM and IL models is
set to be $\tau_{\alpha}=0.01$. For the DPM model, the remaining
hyperparameters are: $\tau_{0}=0.01$, $a_{0}=b_{0}=1$, $f_{0}=r_{\tau}=3$,
$r_{\lambda}=\lambda_{0}=0.03$.

We compared the models in terms of their DIC (to assess overall model fit),
MSE (to assess accuracy in estimating the logROR), as well as the coverage
probability of the resulting 95% posterior credible intervals. We also
generated the caterpillar plots to visualize logROR estimates and heatmaps to
see how well the clustering structure was recovered in the presence of cluster
overlap.

Based on the simulation results in Table 2, the DPM model has dramatically
lower DIC and better MSE compared to the IL-Vague model. The IL-Informative
model performs better than the IL-Vague model (the informative prior helps
with the sparsity present in the data), but its performance is still worse
than the DPM model. All models obtain near nominal coverage for the 95%
posterior credible intervals in the studied cases. We also fit the DPM model
using more informative priors and found the model performance did not change,
suggesting that the model behavior is robust to different hyperparameter
values.

Figure 1: Caterpillar plot from simulation with $\sigma=0.8$. The dashed green
line traces the true $\beta$ values, the dashed red line traces 0, and the 95%
posterior credible intervals are represented by the error bars. Asterisks at
the end of the error bars denote credible intervals outside the range of the
plot.

As shown in Figure 1, the DPM model produces better estimates than the IL
models as evidenced by the proximity of the posterior mean estimates to the
true values, as well as the narrower confidence intervals while maintaining
similar coverage probabilities. For estimates which the DPM model has shrunk
towards 0, the uncertainty in the estimates are properly reflected in the
larger width of the credible intervals. On the other hand, the IL models do
not allow for information sharing between AEs and use the same prior for all
logROR estimates. For the IL-Vague model, we observe more extreme estimates
and much wider confidence intervals than in the DPM model. For the IL-
Informative model, which has a stronger prior to account for the limited
information in the sparse data, while the estimates trace the true values
better, the resulting confidence intervals are noticeably larger than those
from the DPM model. The superior performance in the DPM model is likely due to
the information sharing between similar logROR estimates. Using a mixture of
normal priors while allowing the number of mixture components to be determined
from the data, we are able to induce appropriate shrinkage on the logROR
estimates in a data-adaptive way. The caterpillar plots from the $\sigma=0.5$
and $\sigma=1$ simulations yielded similar conclusions and can be found in
Figures 7 and 8 in Appendix C.

The heatmaps in Figure 2 display the ground truth as well as the learned
clustering structure from the DPM model in the simulation cases. For
$\sigma=0.5$, in which there is no overlap between the clusters, the DPM model
is able to correctly recover the underlying three-cluster structure. In this
case, we attribute the ambiguity in cluster assignment for some of the AEs
primarily to the sparsity in the simulated data, especially for the $\beta$
values coming from the extreme clusters. As $\sigma$ increases and the
between-cluster heterogeneity decreases, the learned clustering structure
becomes less definitive, deviating from the clear block-diagonal structure
seen in the $\sigma=0.5$ simulation. For $\sigma=0.8$, the DPM model is able
to somewhat recover the true cluster structure. We observe more ambiguity in
the cluster assignment from the AEs drawn from $\mu_{k}=0$, however the
extreme clusters are distinct. On the other hand, the DPM model had a much
harder time identifying the cluster structure in the $\sigma=1$ case, in which
we observe significantly more overlap between the center cluster and the
extreme clusters. However, the overall model performance (based on DIC and
MSE) remained better than the IL models, suggesting that the DPM model
allowing information borrowing between logROR parameters has advantages over
the IL models, even in the absence of a clear underlying cluster structure.

(a) Ground Truth $\sigma=0.5$

(b) Ground Truth $\sigma=0.8$

(c) Ground Truth $\sigma=1$

(d) DPM $\sigma=0.5$

(e) DPM $\sigma=0.8$

(f) DPM $\sigma=1$

Figure 2: Ground truth (a-c) and DPM predicted clustering heatmaps (d-f) from
the simulation study. Color indicates the probability that two AEs are
assigned to the same mixture component (i.e. cluster together). The axes
sorted based on the true $\beta$ values in increasing order from left to
right.

##  4 VAERS Data Analysis

In this paper we investigate safety of the mRNA-based COVID-19 vaccines
(Pfizer–BioNTech and Moderna, denoted as “mRNA" vaccines) and the Johnson &
Johnson–Janssen vaccine (denoted “JJ" vaccine). We use influenza vaccines
(denoted as“FLU") as our reference group. FLU includes inactivated and live
attenuated influenza vaccines, labeled as “FLU3" or “FLU4," and “FLUN3" or
“FLUN4" in VAERS.

We fit our proposed model (denoted “DPM" model) and compare results to those
from IL-Informative (i.e, $\tau=0.1$, denoted “IL" model). We did not find
much difference in the results when using the IL-Vague model. We are
interested in the extent to which the information sharing by way of the DPM
prior affects the logROR estimates and model fit. We also applied our proposed
negative control and enrichment methods to identify individual AE signals and
group-level AE signals in the two studies. As described in the Methods
section, an AE is identified as a signal if its NC-adjusted posterior
probability is larger than 90% cutoff and a SOC group $s$ is enriched if the
lower bound of the 95% posterior credible interval of the EOR is greater than
1 and the posterior mean of the EOR greater than 2.

In all studies, we set a maximum of $K=5$ clusters in the DPM model. We use
the same hyperparameter values as in the simulation study for the DPM model.

We ran three independent MCMC chains for each model. For each chain, after a
burn-in of 15,000 iterations we collected one MCMC sample every 5 iterations
until we obtained 1,000 MCMC samples. For model convergence, >99% of the
Gelman-Rubin statistics for the $\beta$ coefficients have $R_{c}$ < 1.2 for
the DPM and IL model MCMC chains.

###  4.1 Data Preprocessing

We conducted our analysis using VAERS data from September 2016 to December
2022. All VAERS records for subjects less than 18 years old were omitted from
the analysis. All VAERS records for COVID-19 vaccines prior to March 16, 2020
were excluded as these occurred prior to the start of COVID-19 vaccine
clinical trials. Finally, all VAERS records indicating administration of both
a COVID-19 and non-COVID-19 vaccine were removed as we would not be able to
differentiate which vaccine was linked to the AE in question.

AEs with less than five total cases were excluded from the analysis. In
addition, we limit our analysis to AEs for which there were both at least one
COVID-19 case and at more than 3 FLU vaccine cases. In the former case, AEs
which have not been reported with a COVID-19 vaccine are not relevant as we
are interested in AEs associated with the COVID-19 vaccine.

We divided subject age into three groups: [18-30), [30-65), and $\geq 65$.
With all categorical variables, we aggregate AE reports by gender, age group,
and vaccine group into binomial count data for computational efficiency.

In the enrichment analysis, we used the SOC level in MedDRA to define groups
of AEs. For AEs belonging to multiple SOC groups—for example “Menstruation
Delayed" belongs to both the “Reproductive system and breast" as well as
“Endocrine" disorder groups—when computing the EOR for one of these groups we
ignore the AEs’ membership in other groups in order to not double-count AEs
and bias the EOR estimates.

###  4.2 Comparing COVID-19 mRNA and FLU vaccines

The final dataset contains 737,301 reports (687,236 mRNA and 50,065 FLU) for
1,089 AEs across 24 SOC groups. Both vaccine groups appear to be predominantly
female and have similar sex distributions; the proportion of Unknown, Female,
and Male subjects were 1.66%, 71.9%, and 26.4% in the FLU group and 1%, 69%,
and 30% in the mRNA group respectively. The FLU group is younger than the mRNA
group, having an age distribution of 40.4%, 49.4%, and 10.2% for the three age
groups compared to 29.0%, 59.4%, and 11.5% for the mRNA group.

The DPM model fit the data better than IL model, obtaining a DIC of 96,242
compared to 96,271 for IL, indicating that appropriate information sharing in
the DPM improves model fit. It uses on average 3.3 mixture components to
estimate the logROR coefficients, where we define a component as “used" if it
has at least 5% of the AEs.

Figure 3: Caterpillar plot in the study for COVID-19 mRNA vaccine, presenting
the 95% posterior credible intervals for the logROR ($\beta$) estimates from
both models. The dashed green line shows values of the posterior means and the
red dashed line corresponds to $\beta=0$.

(a)

(b)

Figure 4: (a) Venn diagram comparing the Signal AE overlap between IL (blue)
and DPM (red) models and (b) Results from SOC group enrichment analysis for
the study of COVID-19 mRNA vaccines.

Figure 3 shows the credible interval estimates for the DPM and IL models. The
IL model credible intervals tend to be wider than those from the DPM model,
which is expected from the information sharing in the DPM prior. In addition,
the more extreme $\beta$ estimates from the IL model tend to be further from 0
than those from the DPM model, especially for the negative $\beta$ estimates
which have a higher incidence of data sparsity and small counts present.

We identified 17 signal AEs for the mRNA vaccines using the DPM model; see
Table 3 in Appendix A for logROR estimates and credible intervals as well as
the NC-adjusted signal probability. Generally the IL and DPM models identified
the same set of signal AEs as shown in the Venn diagram in Figure 4(a). The
differences between the identified signal AEs are that the DPM model
identifies “death," "thrombosis," and “anticoagulant therapy" while IL
identified “myocarditis" and “appendicitis."

The identified signal AEs for pulmonary conditions—acute respiratory failure
and pulmonary embolism—as well as acute myocardial infarction were shown to be
associated with the mRNA COVID-19 vaccines in other studies [30, 31].
Similarly the identified reproductive system signal AEs—irregular and delayed
menstruation as well as menstrual cycle disorders—and vascular system AEs—deep
vein thrombosis and thrombosis—have been linked to the mRNA COVID-19 vaccines
[30, 32, 33]. Acute kidney injury [34], urticaria [35], as well as vaccine
site pruritus and erythema [36, 37, 38] have also been previously reported
with the mRNA COVID-19 vaccines.

Other identified AE signals, including anosmia and ageusia (loss of smell and
taste respectively), death, anticoagulant therapy, and hypoxia (reduced blood
oxygen levels) are associated with COVID-19 infection and/or hospitalization
and the higher reporting rates for these AEs might be related to COVID-19
infection which occurred around the time of vaccination. VAERS takes reports
for AEs that occur after the administration of a vaccine, regardless whether
or not it was caused by the vaccine. Large epidemiological studies are needed
to further interrogate these signals and identify a causal mechanism should
one be present.

Figure 5: Caterpillar plot in the study for the COVID-19 JJ vaccine,
presenting the 95% posterior credible intervals for the logROR ($\beta$)
estimates from both models. The dashed green line shows values of the
posterior means and the red dashed line corresponds to $\beta=0$.

(a)

(b)

Figure 6: (a) Venn diagram comparing the Signal AE overlap between IL (blue)
and DPM (red) models and (b) Results from SOC group enrichment analysis for
the study of the COVID-19 JJ vaccine.

As shown in Figure 4(b), the reproductive system and breast, endocrine,
cardiac, respiratory thoracic and mediastinal, as well as vascular disorder
groups are significantly enriched for the mRNA COVID-19 vaccines. These
results are consistent with the identified signal AEs described above. While
individually AEs may not qualify as signals, if enough AEs within a group are
close to being signals then this can give an ‘overall’ signal for the group of
AEs. The posterior mean EOR values and 95% credible intervals can be found in
Table 5 in Appendix B.

###  4.3 Comparing COVID-19 Johnson & Johnson–Janssen and FLU vaccines

The final dataset contains 106,126 reports (53,901 JJ COVID-19 records and
52,225 FLU) for 1,089 AEs, across 24 SOC groups. Both vaccine groups appear to
be predominantly female, however the JJ group has more balanced female/male
distribution; the proportion of Unknown, Female, Male subjects were 1.61%,
71.6%, and 26.8% in the FLU group and 0.4%, 60.4%, and 39.2% in the mRNA group
respectively. The FLU group is younger than the JJ group, having an age
distribution of 40.2%, 49.6%, and 10.2% for the three age groups,
respectively, compared to 14.5%, 67.7%, and 17.8% for the JJ group.

The DPM model appears to fit the data better than the IL model, obtaining DIC
values of 68,011 and 68,235 respectively.

As in the mRNA vaccine case, Figure 5 shows that IL credible intervals are
wider than those from the DPM model, with larger differences in width for AEs
in the tails of the distribution. The posterior mean logROR estimates appear
to be more extreme for the IL model as well.

We identify 9 signal AEs for the JJ vaccine; see Table 4 in Appendix A for
logROR estimates and credible intervals as well as the NC-adjusted signal
probability. Generally the IL and DPM models identify the same set of signal
AEs, as shown in Figure 6(a), with the differences being the IL model
identifying “thrombectomy" and “menstrual disorder."

Similar to the mRNA vaccine analysis, a host of vascular AEs—pulmonary
embolism, pulmonary thrombosis, deep vein thrombosis, and thrombosis—as well
as AEs related to the reproductive system—menstruation irregular, menstrual
disorder, and dysmenorrhoea. Hypoxia is again identified as a signal AE which,
along with acute respiratory failure, may be related to hospitalization and
the patient disease status rather than the disease itself. Further
epidemiological studies are needed to understand signal AEs and establish any
causal relationships.

We do not observe any significantly enriched SOC groups for the JJ vaccine,
which is likely due to the smaller sample size compared to the mRNA analysis.
The posterior mean EOR values and 95% credible intervals can be found in Table
6 in Appendix B.

##  5 Conclusions

VAERS is an important resource for early identification of vaccine safety
issues. It is not however without limitations, namely small counts for rare
AEs and reporting bias. We propose a Bayesian model using a Dirichlet Process
Mixture (DPM) prior to detect AEs associated with the mRNA-based and Johnson &
Johnson–Janssen COVID-19 vaccines while mitigating the reporting bias. By
applying data-adaptive shrinkage for information sharing, we are able to
improve the estimates for AEs in the presence of data sparsity. Simulation
studies demonstrate the improvements in DIC and MSE from using the proposed
DPM model when compared to a baseline model without information sharing.
Analysis of VAERS data using our proposed model and negative control procedure
identifies previously confirmed AE signals, as well as associated system organ
class groups for the AEs. VAERS collects data on any AE following vaccination,
whether it is related to the vaccine or not. Therefore, we can not draw any
causal conclusion that an AE is caused by a vaccine based on our results. The
vaccine safety signals discovered in our real data analyses need to be
confirmed by larger studies.

Although this paper focuses on the VAERS database, the proposed methods
generally apply to other databases which rely on passive reporting, such as
FDA Adverse Events Reporting System (FAERS) and the Adverse Drug Reactions
(ADR) database for conducting post-marketing drug safety surveillance.

Supplementary Material

  * •

Supplementary Document A: JAGS Code for DPM Model

  * •

Supplementary Document B: List of VAERS NC AEs

Funding

Research reported in this publication was supported by the National Institute
Of Allergy And Infectious Diseases of the National Institutes of Health under
Award Number R01AI158543. The content is solely the responsibility of the
authors and does not necessarily represent the official views of the National
Institutes of Health.

## References

  * [1] Edouard Mathieu, Hannah Ritchie, Esteban Ortiz-Ospina, Max Roser, Joe Hasell, Cameron Appel, Charlie Giattino, and Lucas Rodés-Guirao.  A global database of covid-19 vaccinations.  Nature human behaviour, 5(7):947–953, 2021.
  * [2] Gianmarco Troiano and Alessandra Nardi.  Vaccine hesitancy in the era of covid-19.  Public Health, 2021.
  * [3] Tom T Shimabukuro, Michael Nguyen, David Martin, and Frank DeStefano.  Safety monitoring in the vaccine adverse event reporting system (vaers).  Vaccine, 33(36):4398–4405, 2015.
  * [4] Tom T Shimabukuro, Matthew Cole, and John R Su.  Reports of anaphylaxis after receipt of mrna covid-19 vaccines in the us—december 14, 2020-january 18, 2021.  Jama, 325(11):1101–1102, 2021.
  * [5] Matthew E Oster, David K Shay, John R Su, Julianne Gee, C Buddy Creech, Karen R Broder, Kathryn Edwards, Jonathan H Soslow, Jeffrey M Dendy, Elizabeth Schlaudecker, et al.  Myocarditis cases reported after mrna-based covid-19 vaccination in the us from december 2020 to august 2021.  JAMA, 327(4):331–340, 2022.
  * [6] Emily Jane Woo, Adamma Mba-Jonas, Rositsa B Dimova, Meghna Alimchandani, Craig E Zinderman, and Narayan Nair.  Association of receipt of the ad26. cov2. s covid-19 vaccine with presumptive guillain-barré syndrome, february-july 2021.  JAMA, 2021.
  * [7] Lili Zhao, Sunghun Lee, Rongxia Li, Edison Ong, Yongqun He, and Gary Freed.  Improvement in the analysis of vaccine adverse event reporting system database.  Statistics in biopharmaceutical research, 12(3):303–310, 2020.
  * [8] William DuMouchel.  Bayesian data mining in large frequency tables, with an application to the fda spontaneous reporting system.  The American Statistician, 53(3):177–190, 1999.
  * [9] S. J. Evans, P. C. Waller, and S. Davis.  Use of proportional reporting ratios (PRRs) for signal generation from spontaneous adverse drug reaction reports.  Pharmacoepidemiology and Drug Safety, 10:483–486, 2001.
  * [10] Lan Huang, Jyoti Zalkikar, and Ram C Tiwari.  A likelihood ratio test based method for signal detection with application to fda’s drug safety data.  Journal of the American Statistical Association, 106(496):1230–1241, 2011.
  * [11] Lan Huang, Ted Guo, Jyoti N Zalkikar, and Ram C Tiwari.  A review of statistical methods for safety surveillance.  Therapeutic Innovation & Regulatory Science, 48(1):98–108, 2014\.
  * [12] E. P. van Puijenbroek, A. Bate, H. G. Leufkens, M. Lindquist, R. Orre R, and A. C. Egberts.  A comparison of measures of disproportionality for signal detection in spontaneous reporting systems for adverse drug reactions.  Pharmacoepidemiology and Drug Safety, 11:3–10, 2002.
  * [13] A. Bate, M. Lindquist, I. R. Edwards, S. Olsson, R. Orre, A. Lansner, and R. M. De Freitas.  A Bayesian neural network method for adverse drug reaction signal generation.  European Journal of Clinical Pharmacology, 54:315–321, 1998.
  * [14] R. Orre A. Lansner, A. Bate, I. R. Edwards M. Lindquist, S. Olsson, R. Orre, A. Lansner, and R. M. De Freitas.  Bayesian neural networks with confidence estimations applied to data mining.  Comput Stat Data Anal, 34:473–493, 2000.
  * [15] G. N. Nóren, A. Bate, R. Orre, and I. R. Edwards.  Extending the methods used to screen the WHO drug safety database towards analysis of complex associations and improved accuracy for rare events.  BCPNN, 25:3740–3757, 2006.
  * [16] W. DuMouchel and D. Pregibon.  Empirical bayes screening for multi-item associations.  Proceedings of the Seventh ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, San Francisco,CA, page 67–76, 2001\.
  * [17] A. Szarfman, S. G. Machado, and R. T. O’Neill.  Use of screening algorithms and computer systems to efficiently signal higher-than-expected combinations of drugs and events in the us FDA’s spontaneous reports database.  PLoS ONE, 25:381–392, 2002.
  * [18] M. Kulldorff, R. L. Davis, M. Kolczak, E. Lewis, T. Lieu, and R. Platt.  A maximized sequential probability ratio test for drug and vaccine safety surveillance.  Sequential Analysis, 30:58–78, 2011.
  * [19] R. L. Davis, M. Kolczak, E. Lewis, J. Nordin, M. Goodman, D. K. Shay, R. Platt, S. Black, H. Shinefield, and R. T Chen.  Active surveillance of vaccine safety: a system to detect early signs of adverse events.  Epidemiology, 16:336–341, 2005.
  * [20] L. Li and M. Kulldorff.  A conditional maximized sequential probability ratio test for pharmacovigilance.  Statistics in Medicine, 29:284–295, 2009.
  * [21] Rongxia Li, Brock Stewart, Eric Weintraub, , and Michael M. McNeil.  Continuous sequential boundaries for vaccine safety surveillance.  Statistics in Medicine, 33:3387–3397, 2014.
  * [22] M Kulldorff, I. Dashevsky, T. R. Avery, K. A. Chan, R. L. Davis, D. Graham, R. Platt, S. E. Andrade, D. Boudreau, M. J. Gunter, L. J. Herrinton, P. Pawloski, M. A. Raebel, D. Roblin, and J. S. Brown.  Drug safety data mining with a tree-based scan statistic.  Pharmacoepidemiology and Drug Safety, 22:517–523, 2013.
  * [23] David I Ohlssen, Linda D Sharples, and David J Spiegelhalter.  Flexible random-effects models using bayesian semi-parametric models: applications to institutional comparisons.  Statistics in medicine, 26(9):2088–2112, 2007.
  * [24] Xu Shi, Wang Miao, and Eric Tchetgen.  A selective review of negative control methods in epidemiology.  Current Epidemiology Reports, 7, 10 2020.
  * [25] Martijn J Schuemie, Patrick B Ryan, William DuMouchel, Marc A Suchard, and David Madigan.  Interpreting observational studies: why empirical calibration is needed to correct p-values.  Statistics in medicine, 33(2):209–218, 2014.
  * [26] E. G. Brown, L. Wood, and S. Wood.  The medical dictionary for regulatory activities (meddra) drug safety.  The American Statistician, 20:109–117, 1999.
  * [27] Jayaram Sethuraman.  A constructive definition of dirichlet priors.  Statistica sinica, pages 639–650, 1994.
  * [28] Hemant Ishwaran and Mahmoud Zarepour.  Markov chain monte carlo in approximate dirichlet and beta two-parameter process hierarchical models.  Biometrika, 87(2):371–390, 2000.
  * [29] Andrew Gelman.  Prior distributions for variance parameters in hierarchical models (comment on article by browne and draper).  Bayesian Analysis, 1(3):515 – 534, 2006.
  * [30] Yusuke Miyazato, Kei Yamamoto, Gen Yamada, Shuji Kubota, Masahiro Ishikane, Masaya Sugiyama, Mikako Ueno, Akihiro Matsunaga, Tohru Miyoshi-Akiyama, Yukihito Ishizaka, et al.  Multisystem inflammatory syndrome in adult after first dose of mrna vaccine.  Emerging infectious diseases, 28(4):870, 2022.
  * [31] Marie Joelle Jabagi, Jérémie Botton, Marion Bertrand, Alain Weill, Paddy Farrington, Mahmoud Zureik, and Rosemary Dray-Spira.  Myocardial infarction, stroke, and pulmonary embolism after bnt162b2 mrna covid-19 vaccine in people aged 75 years or older.  Jama, 327(1):80–82, 2022.
  * [32] Alison Edelman, Emily R Boniface, Eleonora Benhar, Leo Han, Kristen A Matteson, Carlotta Favaro, Jack T Pearson, and Blair G Darney.  Association between menstrual cycle length and coronavirus disease 2019 (covid-19) vaccination: a us cohort.  Obstetrics & Gynecology, 139(4):481–489, 2022.
  * [33] T Pietri, J Micallef, M Gervoise-Boyer, and P Boyer.  P-308 spontaneous reports of menstrual cycle disorders after mrna covid-19 vaccine.  Human Reproduction, 37(Supplement_1):deac107–294, 2022.
  * [34] Vivette D D’Agati, Satoru Kudose, Andrew S Bomback, Ananea Adamidis, and Albert Tartini.  Minimal change disease and acute kidney injury following the pfizer-biontech covid-19 vaccine.  Kidney international, 100(2):461–463, 2021.
  * [35] Eli Magen, Avi Yakov, Ilan Green, Ariel Israel, Shlomo Vinker, and Eugene Merzon.  Chronic spontaneous urticaria after bnt162b2 mrna (pfizer-biontech) vaccination against sars-cov-2.  In Allergy and Asthma Proceedings, volume 43, page 30. OceanSide Publications, 2022.
  * [36] C Doyle and Richard E Watchorn.  Persistent pruritus following covid-19 vaccination.  Irish Journal of Medical Science (1971-), pages 1–2, 2022.
  * [37] BM Cyrenne, F Al-Mohammedi, JG DeKoven, and R Alhusayen.  Pityriasis rosea-like eruptions following vaccination with bnt162b2 mrna covid-19 vaccine.  Journal of the European Academy of Dermatology and Venereology, 2021\.
  * [38] Luca Borg, Liam Mercieca, Dillon Mintoff, Daniel Micallef, David Pisani, Alexandra Betts, and Lawrence Scerri.  Pfizer-biontech sars-cov-2 mrna vaccine-associated erythema multiforme.  Journal of the European Academy of Dermatology and Venereology, 36(1):e22–e24, 2022.

##  Appendix A Signal AE Tables

Adverse Event Name |  $\beta$ Post. Mean |  $\beta$ 95% Credible Interval | Signal Prob.
---|---|---|---
Menstruation irregular | 3.47 | (2.72, 4.41) | 1.00
Menstrual disorder | 3.41 | (2.78, 4.17) | 1.00
Vaccination site pruritus | 3.11 | (2.67, 3.53) | 1.00
Deep vein thrombosis | 2.82 | (2.26, 3.46) | 1.00
Acute myocardial infarction | 2.78 | (2.12, 3.52) | 1.00
Hypoxia | 2.69 | (2.34, 3.06) | 1.00
Acute respiratory failure | 2.61 | (2.25, 3.00) | 1.00
Anosmia | 2.57 | (2.13, 2.96) | 1.00
Vaccination site erythema | 2.31 | (2.07, 2.58) | 1.00
Pulmonary embolism | 2.69 | (2.20, 3.18) | 1.00
Ageusia | 2.33 | (1.96, 2.76) | 0.99
Anticoagulant therapy | 2.11 | (1.79, 2.44) | 0.96
Acute kidney injury | 2.22 | (1.79, 2.64) | 0.96
Mechanical urticaria | 2.43 | (1.69, 3.91) | 0.96
Death | 1.98 | (1.82, 2.14) | 0.94
Menstruation delayed | 2.23 | (1.53, 3.07) | 0.91
Thrombosis | 2.00 | (1.65, 2.42) | 0.91
Table 3: Adverse events with at least a 90% signal probability for the mRNA COVID-19 vaccines when compared to the FLU vaccine. Posterior mean of the logROR and the 95% credible interval are also presented. Adverse Event Name |  $\beta$ Post. Mean |  $\beta$ 95% Credible Interval | Signal Prob.
---|---|---|---
Vaccine breakthrough infection | 3.84 | (3.41, 4.30) | 1.00
Deep vein thrombosis | 3.49 | (3.08, 3.97) | 1.00
Pulmonary embolism | 3.05 | (2.70, 3.43) | 1.00
Menstruation irregular | 3.16 | (2.64, 3.82) | 1.00
Thrombosis | 2.89 | (2.59, 3.18) | 0.99
Hypoxia | 2.84 | (2.52, 3.18) | 0.99
Acute respiratory failure | 2.79 | (2.47, 3.14) | 0.97
Dysmenorrhoea | 3.03 | (2.34, 3.75) | 0.97
Pulmonary thrombosis | 2.99 | (2.28, 3.87) | 0.94
Table 4: Adverse events with at least a 90% signal probability for the Johnson
and Jonson COVID-19 vaccine when compared to the FLU vaccine. Posterior mean
of the logROR and the 95% credible interval are also presented.

##  Appendix B SOC Group Enrichment Tables

Name | EOR | 95% Credible Interval
---|---|---
Reproductive system and breast disorders | 9.66 | (4.77, 18.62)
Endocrine disorders | 5.01 | (1.93, 10.47)
Cardiac disorders | 3.94 | (1.02, 7.87)
Respiratory thoracic and mediastinal disorders | 2.64 | (1.59, 4.15)
Vascular disorders | 2.15 | (1.17, 3.59)
Metabolism and nutrition disorders | 1.88 | (0.70, 4.20)
Surgical and medical procedures | 1.83 | (0.78, 3.37)
Renal and urinary disorders | 1.56 | (0.00, 3.70)
Injury poisoning and procedural complications | 1.37 | (0.82, 2.04)
General disorders and administration site conditions | 0.97 | (0.57, 1.49)
Hepatobiliary disorders | 0.96 | (0.00, 13.55)
Skin and subcutaneous tissue disorders | 0.71 | (0.38, 1.30)
Pregnancy puerperium and perinatal conditions | 0.58 | (0.00, 1.73)
Gastrointestinal disorders | 0.55 | (0.00, 1.18)
Infections and infestations | 0.47 | (0.00, 1.31)
Immune system disorders | 0.47 | (0.00, 1.15)
Nervous system disorders | 0.35 | (0.16, 0.63)
Psychiatric disorders | 0.18 | (0.00, 0.68)
Social circumstances | 0.09 | (0.00, 1.78)
Blood and lymphatic system disorders | 0.04 | (0.00, 0.00)
Eye disorders | 0.02 | (0.00, 0.32)
Musculoskeletal and connective tissue disorders | 0.01 | (0.00, 0.18)
Ear and labyrinth disorders | 0.01 | (0.00, 0.00)
Neoplasms benign malignant and unspecified incl cysts and polyps | 0.00 | (0.00, 0.00)
Table 5: SOC group enrichment OR (EOR) and 95% credible intervals for COVID-19 mRNA vs FLU vaccine analysis. Name | EOR | 95% Credible Interval
---|---|---
Reproductive system and breast disorders | 9.27 | (0.00, 31.08)
Vascular disorders | 7.58 | (0.00, 17.64)
Endocrine disorders | 5.99 | (0.00, 19.91)
Infections and infestations | 4.43 | (0.00, 14.77)
Respiratory thoracic and mediastinal disorders | 2.24 | (0.00, 6.28)
General disorders and administration site conditions | 1.71 | (0.42, 5.88)
Injury poisoning and procedural complications | 1.64 | (0.00, 6.06)
Surgical and medical procedures | 1.35 | (0.00, 7.74)
Metabolism and nutrition disorders | 0.86 | (0.00, 4.78)
Cardiac disorders | 0.43 | (0.00, 3.37)
Renal and urinary disorders | 0.20 | (0.00, 3.11)
Nervous system disorders | 0.12 | (0.00, 0.81)
Gastrointestinal disorders | 0.09 | (0.00, 1.06)
Hepatobiliary disorders | 0.07 | (0.00, 0.00)
Blood and lymphatic system disorders | 0.03 | (0.00, 0.00)
Neoplasms benign malignant and unspecified incl cysts and polyps | 0.03 | (0.00, 0.00)
Psychiatric disorders | 0.01 | (0.00, 0.00)
Eye disorders | 0.01 | (0.00, 0.00)
Social circumstances | 0.01 | (0.00, 0.00)
Musculoskeletal and connective tissue disorders | 0.00 | (0.00, 0.00)
Immune system disorders | 0.00 | (0.00, 0.00)
Skin and subcutaneous tissue disorders | 0.00 | (0.00, 0.00)
Ear and labyrinth disorders | 0.00 | (0.00, 0.00)
Pregnancy puerperium and perinatal conditions | 0.00 | (0.00, 0.00)
Table 6: SOC group enrichment OR (EOR) and 95% credible intervals for COVID-19
JJ vs FLU vaccine analysis.

##  Appendix C Simulation Caterpillar Plots

Figure 7: Caterpillar plot from simulation with $\sigma=0.5$. The dashed green
line traces the true $\beta$ values, the dashed red line traces 0, and the 95%
posterior credible intervals are represented by the error bars. Asterisks at
the end of the error bars denote credible intervals outside the range of the
plot. Figure 8: Caterpillar plot from simulation with $\sigma=1$. The dashed
green line traces the true $\beta$ values, the dashed red line traces 0, and
the 95% posterior credible intervals are represented by the error bars.
Asterisks at the end of the error bars denote credible intervals outside the
range of the plot.