summaries: This paper describes a method for the causal analysis of long streams of event data that can be modeled as multivariate point processes (MPPs). The paper defines the average treatment effect (ATE) between pairs of events in MPPs that occur within a window, and derives how to estimate the ATE as the expected difference of average outcome rates.  When covariates differ or assumptions are violated, this estimate may be biased, but adjusting for the propensity score yields an unbiased estimator, which the paper does through inverse probability of treatment weighting (IPTW) (including a stabilized version).  Experiments on data from known, synthetic MPPs (proximal graphical event models and Hawkes processes) show that the proposed IPTW estimation method more accurately recovers the ATE than baseline conditional intensity scores.  Similarly, experiments on a diabetes dataset show that stable IPTW recovers more true causal pairs at 2 out of 3 levels of recall compared to the conditional intensity scores.  These results extend to using the cumulative rather than the conditional intensity rate.

summaries: This paper proposes a model for estimating the average causal effect in point processes where the estimand of interest is the effect of an event occurring within a window on the rate parameter of a point process. The authors use the formulation commonly used for inverse propensity scores of generalized treatments (e.g., from Robins, Imai). The authors describe the properties of the propensity score as a balancing weight by leveraging known results and advocate for the use of stabilization of inverse propensity scores to reduce variance as in Robins(2000). The weights are then used within an outcome model which is a point process model. Empirical results show promising performance compared to a range of alternatives.

summaries: This paper proposes a causal inference framework for a multivariate point process model. The key idea is to form a series of binary processes Z_t and X_t using a moving window [t-w,t] to indicate whether processes Z and X have at least one event in the window [t-w,t].  An average treatment event (ATE) is defined and propensity scores are proposed for unbiased estimation of ATE. The authors use synthetic experiments and real-world datasets to demonstrate how the proposed method work compared to baseline approaches.

summaries: The paper provides tools to estimate average treatment effect (ATE) between event pairs. The authors first define the quantity for multivariate point processes together with the necessary concepts from Neyman–Rubin causality, then they derive a propensity score in this setting. Finally, they give a method for estimating the ATE and evaluate the method on synthetic and real world datasets.

abstractText: Causal inference and discovery from observational data has been extensively studied across multiple fields. However, most prior work has focused on independent and identically distributed (i.i.d.) data. In this paper, we propose a formalization for causal inference between pairs of event variables in multivariate recurrent event streams by extending Rubin’s framework for the average treatment effect (ATE) and propensity scores to multivariate point processes. Analogous to a joint probability distribution representing i.i.d. data, a multivariate point process represents data involving asynchronous and irregularly spaced occurrences of various types of events over a common timeline. We theoretically justify our point process causal framework and show how to obtain unbiased estimates of the proposed measure. We conduct an experimental investigation using synthetic and real-world event datasets, where our proposed causal inference framework is shown to exhibit superior performance against a set of baseline pairwise causal association scores.

authors:
Odd O Aalen
Mats J Stensrud
Vanessa Didelez
Rhian Daniel
Kjetil Røysland
Susanne Strohmaier

title: Time-dependent mediators in survival analysis: Modeling direct and indirect effects 1DISTRIBUTION STATEMENT A. Approved for public release: distribution unlimited 10 with the additive hazards model

venue: Biometrical Journal,

title: Causal Modeling and Prediction over Event Streams

venue: PhD thesis, University of Vermont,

authors:
V. Arya
K. Shanmugam
P. Aggarwal
Q. Wang
P. Mohapatra
S. Nagar

title: Evaluation of causal inference techniques for AIOps

venue: In 8th ACM IKDD CODS and 26th COMAD,

authors:
Serge Assaad
Shuxi Zeng
Chenyang Tao
Shounak Datta
Nikhil Mehta
Ricardo Henao
Fan Li
Lawrence Carin Duke

title: Counterfactual representation learning with balancing weights

venue: In International Conference on Artificial Intelligence and Statistics,

authors:
Susan Athey
Guido W Imbens
Stefan Wager

title: Approximate residual balancing: debiased inference of average treatment effects in high dimensions

venue: Journal of the Royal Statistical Society: Series B (Statistical Methodology),

authors:
E. Bacry
M. Bompaire
S. Gaïffas
S. Poulsen

title: tick: A Python library for statistical learning, with a particular emphasis on time-dependent modeling

authors:
D. Bhattacharjya
D. Subramanian
T. Gao

title: Proximal graphical event models

venue: In Advances in Neural Information Processing Systems (NeurIPS),

authors:
Debarun Bhattacharjya
Tian Gao
Nicholas Mattei
Dharmashankar Subramanian

title: Causeeffect association between event pairs in event datasets

venue: In International Joint Conference on Artificial Intelligence (IJCAI),

authors:
D.R. Cox
P.A.W. Lewis

title: Multivariate point processes

venue: In Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 3: Probability Theory,

authors:
Daryl J Daley
D Vere Jones

title: An Introduction to the Theory of Point Processes

venue: Elementary Theory of Point Processes. Springer,

title: Graphical models for marked point processes based on local independence

venue: Journal of the Royal Statistical Society, Series B,

authors:
Michael Eichler
Rainer Dahlhaus
Johannes Dueck

title: Graphical modeling for multivariate Hawkes processes with nonparametric link functions

venue: Journal of Time Series Analysis,

authors:
Christian Fong
Chad Hazlett
Kosuke Imai

title: Covariate balancing propensity score for a continuous treatment: Application to the efficacy of political advertisements

venue: The Annals of Applied Statistics,

authors:
A. Gunawardana
C. Meek

title: Universal models of multivariate temporal point processes

venue: In Proc. of the 19th AISTATS,

authors:
A. Gunawardana
C. Meek
P. Xu

title: A model for temporal dependencies in event streams

authors:
Ruocheng Guo
Jundong Li
Huan Liu

title: Initiator: Noise-contrastive estimation for marked temporal point process

venue: In Proceedings of the International Joint Conference in Artificial Intelligence (IJCAI),

authors:
Guido W Imbens
Donald B Rubin

title: Causal inference in statistics, social, and biomedical sciences

authors:
Fredrik Johansson
Uri Shalit
David Sontag

title: Learning representations for counterfactual inference

venue: In Proceedings of the International Conference on Machine Learning (ICML),

authors:
Quyu Kong
Marian-Andrei Rizoiu
Lexing Xie

title: Exploiting uncertainty in popularity prediction of information diffusion cascades using self-exciting point processes

venue: arXiv preprint arXiv:2001.11132,

authors:
Fan Li
Kari Lock Morgan
Alan M Zaslavsky

title: Balancing covariates via propensity score weighting

venue: Journal of the American Statistical Association,

title: Statistical modeling of causal effects in continuous time

venue: The Annals of Statistics,

authors:
Huzhang Mao
Liang Li
Wei Yang
Yu Shen

title: On the propensity score weighting analysis with survival outcome: Estimands, estimation, and inference

venue: Statistics in Medicine,

authors:
Hongyuan Mei
Jason M Eisner

title: The neural Hawkes process: A neurally self-modulating multivariate point process

title: Optimal dynamic treatment regimes

title: Crisis early warning and decision support: Contemporary approaches and thoughts on future research

venue: International Studies Review,

authors:
Y. Okada
K.I. Fukui
K. Moriyama
M. Numao

title: Causal inference with time and space proximity under uncertainty

venue: In Proceedings of the Pacific-Asia Conference on Knowledge Discovery and Data Mining,

authors:
A.P. Parikh
A. Gunawardana
C. Meek

title: Conjoint modeling of temporal dependencies in event streams

venue: In Proceedings of Uncertainty in Artificial Intelligence Workshop on Bayesian Modeling Applications,

title: Causality: Models, Reasoning and Inference

authors:
S. Rajaram
T. Graepel
R. Herbrich

title: Poisson-networks: A model for structured point processes

venue: In Proceedings of the Conference on Artificial Intelligence and Statistics (AISTATS),

title: Temporal point processes: the conditional intensity function

venue: Lecture Notes,

authors:
Marian-Andrei Rizoiu
Swapnil Mishra
Quyu Kong
Mark Carman
Lexing Xie

title: Sir- Hawkes: linking epidemic models and Hawkes processes to model diffusions in finite populations

venue: In Proceedings of the World Wide Web (WWW) Conference,

authors:
James M Robins
Miguel Angel Hernan
Babette Brumback

title: Marginal structural models and causal inference in epidemiology

authors:
Paul R Rosenbaum
Donald B Rubin

title: The central role of the propensity score in observational studies for causal effects

title: Counterfactual analyses with graphical models based on local independence

authors:
Pål C Ryalen
Mats J Stensrud
Kjetil Røysland

title: Transforming cumulative hazard

venue: estimates. Biometrika,

authors:
Pål Christie Ryalen
Mats Julius Stensrud
Sophie Fosså
Kjetil Røysland

title: Causal inference in continuous time: an example on prostate cancer

venue: therapy. Biostatistics,

authors:
Peter Schulam
Suchi Saria

title: Reliable decision support using counterfactual models

title: Discretizing logged interaction data biases learning for decisionmaking

venue: arXiv preprint arXiv:1810.03025,

authors:
Rajat Sen
Hsiang-Fu Yu
Inderjit S. Dhillon

title: Think globally, act locally: A deep neural network approach to high-dimensional time series forecasting

venue: In Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems

authors:
Oleg Sofrygin
Zheng Zhu
Julie A Schmittdiel
Alyce S Adams
Richard W Grant
Mark J van der Laan
Romain Neugebauer

title: Targeted learning with daily ehr data

venue: Statistics in medicine,

authors:
J.C. Weiss
D. Page

title: Forest-based point process for event prediction from electronic health records

venue: In Machine Learning and Knowledge Discovery in Databases,

authors:
Shuai Xiao
Junchi Yan
Xiaokang Yang
Hongyuan Zha
Stephen M Chu

title: Modeling the intensity function of point process via recurrent neural networks

venue: In Proceedings of the Conference on Artificial Intelligence (AAAI),

authors:
Shu Yang
Karen Pieper
Frank Cools

title: Semiparametric estimation of structural failure time models in continuous-time processes

authors:
Mingyuan Zhang
Marshall M Joffe
Dylan S Small

title: Causal inference for continuous-time processes when covariates are observed only at discrete times

venue: Annals of Statistics,

authors:
Qiang Zhang
Aldo Lipani
Omer Kirnap
Emine Yilmaz

title: Self-attentive hawkes processes

venue: arXiv preprint arXiv:1907.07561,

heading: 1 Introduction

text:
It is widely known that the gold standard for effective causal inference is through the use of interventional data such as randomized controlled trials, for measuring the impact of some treatment on an outcome of interest [30]. However, intervening in the system can often be impractical or even impossible, in which case one must conduct causal analysis using observational data alone.
Observational data in the form of multivariate event streams is common in many domains including health, finance, and retail. Analogous to how an observation in an i.i.d dataset can be viewed as a sample from a joint distribution over a set of random variables, a multivariate event stream can be viewed as a sample from a multivariate point process over a set of event labels [10], where the instantaneous rate of any event label depends on the historical occurrences of some subset of the event labels. Modeling, fitting, and predicting future occurrences given a multivariate event stream is an active area in statistics [9, 11, 40], data mining [28, 34], and machine learning [32, 16, 29, 44, 7]. Recent work in this area leverages advances in deep learning [45, 25, 17, 48, 21].
Related work in survival analysis focuses only on pairs of events with continuous covariates [23], involving short event streams with both continuous and discrete variables, and most notably [47, 46] only a single occurrence of the outcome [24]. While a general theory for point process counterfactual inference has been developed [37, 1], practical algorithms focus on continuous covariates and hazard models [39]. Moreover, work in dynamic treatment [26, 18] is for discrete time, motivating work for continuous-time data, which aims to reduce bias or high variance in discretizing time [38, 41, 43]. In
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
contrast, here we focus on the setting where data is of the form of long event streams with multiple occurrences of various types of events, including both the treatment event and outcome event.
Causal inference entails drawing a conclusion about a causal relationship between potential causes and effects; in this work, we focus on causal inference between pairs of event labels observed in event stream data, which differs from standard Granger causal graph learning [12, 49]. Specifically, we pose the following causal inference problem for multivariate point processes: how can one meaningfully measure the causal relationship between a cause event label z and an effect event label y? Such a causal measure would explain whether the cause label z amplifies, inhibits, or has no impact on the effect label y, while taking into account the potential effects of all other event labels x.
This problem deviates from the typical causal inference setting in several ways. First, most methods assume i.i.d. observations; in our setting, events may be correlated across time and therefore it is necessary to make certain independence assumptions in order to identify the effect of interest. Second, most causal inference settings are concerned with the expected value of a single observable outcome such as mortality. Our setting has repeated occurrences of an outcome over time and therefore involves event frequency (or rate), which existing pairwise causal inference models are unable to handle and therefore necessitates new methods. Note that since event label z’s historical occurrences could have a complex dynamic effect on another event y’s occurrence at any time in a multivariate point process, there are a number of questions to address around the fundamental problem of even defining the notions of treatment, outcome and covariates, before one can study the causal effect of z on y.
There are several real-world applications that motivate our work on causal inference with event stream data. For instance, one may be interested in understanding causal relations between events relevant to a diabetic patient, such as meals, insulin intake, exercise and blood glucose level changes [2]. Detecting issues in microservice or other web applications requires detecting which faults or abnormal states cause other downstream faults from system logs [3]. Example causal pairs from other domains are: food shortage and protests, which are recurring event labels in socio-economic datasets [27], and earning call events and stock price jumps, which are recurring event labels in financial datasets [42]. All these examples have an observable outcome that occurs repeatedly, necessitating a new approach.
Contributions We make several major contributions: 1) we propose a formalization of causal inference between a pair of event variables in a multivariate point process, 2) we define an average treatment effect (ATE) for multivariate point processes as the effect of having observed treatment in a proximal window, 3) we derive and obtain the equivalence of propensity scores and balancing scores for ATE in multivariate point processes, 4) we propose an inverse propensity score re-weighting procedure for better ATE estimation to adjust for the impact of other covariates, and 5) we conduct a detailed empirical evaluation on both synthetic and real datasets with baselines to show the superior performance of our proposed methods in pairwise causal inference in multivariate point processes.

heading: 2 Background

text:
A multivariate event stream (or event dataset) is a sequence of events, D = {Di}Ni=1, where each event Di includes an event label (or synonymously type or variable) and a time-stamp, i.e. Di = (xi, ti). xi belongs to the label set L, whose cardinality is M , i.e. |L| = M , and ti is the time of occurrence, ti 2 R+. We assume a strictly temporally ordered dataset, ti < tj for i < j, initial time t0 = 0 and end time tN+1 = T . Let z, y refer to an arbitrary pair of event labels in L.
Multivariate event streams can be regarded as samples from a multivariate point process (MPP), where each label has a counting process [10]. It uses conditional intensity functions y(t|Ht) > 0 that denote the instantaneous rate of type y given the history up to t, i.e. Ht = {(xi, ti) : ti < t}. In a multivariate point process, the probability of observing y as the next event at time t is:
pt(y|Htn) = y(t|Ht) exp(
Z t
tn
X
x2L
x(⌧ |H⌧ )d⌧) (1)
where tn is the most recent event occurrence time before t. It can be shown that y(t|Ht)dt = E[Ny([t, t+ dt])|Ht] [10, 33], where Ny(A) is the number of y occurrences in time interval A.
Prior work has proposed the notion of process independence among event labels’ counting processes [11]. The basic idea is that the intensity of one type of event does not depend on certain past events once we know about specific other past events. This is an asymmetric concept, similar to Granger causality. Informally, for sets of labels x, y and z s.t. y \ z = ;, x is process independent of y given
z when all labels in x have conditional intensities that do not functionally depend at any time on the history of labels in y, given the history of labels in z. Process independence is captured in graphical event models [11, 15], where the conditional intensity for any event label at any time t depends only on historical occurrences of its event label parents in the underlying graph.
A minimal graphical representation can be used to define direct causes in multivariate point processes, analogous to those in causal networks. Here we provide a simplified definition of direct causes [11]: Definition 1. Event label z is a direct cause of label y if z belongs to the minimal set of nodes u s.t. y is process independent of all other labels given u.
A minimal set entails that none of its subsets satisfies the property. We highlight that in this work, we are interested in identifying pairwise causal relationships without learning a full graphical model. This is because it can often be computationally intensive as well as unnecessary to learn a full model when one is merely interested in pairwise causality.

heading: 3 Causal Inference in Point Processes

text:
The Neyman-Rubin potential outcomes causal inference framework estimates the treatment effects between a treated variable and an effect for i.i.d. data [36, 20]. We extend this framework, including average treatment effects (ATE), to study how an event label z affects event label y in a MPP. In this class of models, treatment is denoted as z 2 {0, 1}, where 0 is the control and 1 the treatment. The potential outcome yz for each z is the outcome if treatment z is applied. The main difficulty in causal inference stems from the fact that only outcomes from the administered treatments are observed and never any of the other outcomes. Hence it may be viewed as a missing data problem.
We estimate the treatment effects between a treatment variable associated with historical occurrences of z and an outcome variable (or response) associated with y’s occurrence, under existing assumptions in the literature. Covariate variables x involve historical occurrences of labels other than z, i.e. x = L \ z. We assume all variables are observed, hence strong ignorability is satisfied (which we will discuss further in Section 3.2). Next we define treatment, outcome, and covariate variables, and then derive the propensity score for ATE computation.

heading: 3.1 Defining Treatment, Outcome and Covariates

text:
There are many possible ways to define the treatment, outcome, and covariates in a multivariate point process. We begin with the following general formulation. We use Zt to represent the treatment variable at time t, distinguishing it from the event label z. Let Hzt be the history of occurrences of event z before t, and Kyt be the future of event y occurrences at t. Definition 2. General Formulation for Causal Inference: for a pair of event labels (z, y), the treatment variable Zt at time t is a function of z’s historical occurrences, i.e. Zt = fz(Hzt ). The outcome variable is a function of y’s future occurrences at t, i.e. Yt = fy(K y t ). The covariates Xt are a function of historical occurrences of event labels other than z, i.e. Xt = f z(H z t ).
Treatment and Covariate Models We need to specify the functions in the general formulation of Definition 2, for instance, we need to summarize historical occurrences Hzt . Due to potential multiple (re)-occurrence of treatment and outcome, the representation of history needs careful modeling based on specific assumptions at each time t. We take a simplified yet practical view of treatment, making the proximal assumption that the most recent history is sufficient for causation. We do this by assuming only a recent time window [t w, t) for some window w could directly have a causal impact on the outcome at time t; such an assumption has been used previously in point process models [7] and provides an efficient and interpretable definition for treatment. It fits many real-world situations where the influence comes primarily from recent history and repeated occurrences have little additional impact on the intensity rate. It also has a strong statistical justification: matching long treatment history could result in infinite sample requirements. We assume covariates are also defined by whether they have appeared at least once in some recent time window.
Outcome Model We wish to measure a statistical quantity fy(Kyt ) associated with the outcome in a multivariate point process, based on future occurrences of y, Kyt , from time t. Many such quantities fy(Kyt ) can be used. It is natural to consider either the instantaneous expected number of occurrences of the outcome y, E[Ny([t, t+ dt])|Ht], in place of Kyt , or cumulative expected number
of occurrences of outcome in a future duration wf , E[Ny([t, t+ wf ])|Ht]. Moreover, the expected count E[Ny([t, t+ dt])|Ht] equals y(t)dt by definition. Note that at each time, the history can be very different with the proximal assumption, hence the rate at each time in Kyt can also be different. The cumulative counts can be similarly written as ⇤wfy (t) = R t+wf t y(t)dt for some future duration wf , which can be the total time horizon or infinity, but practically it is more meaningful for a specific duration suitable for the application.
We will focus our analysis with y(t) as the outcome in the following, which could be directly extended to outcome ⇤wfy (t). We will compare with both definitions of outcomes in experiments.
Hence, our specific realization of the general formulation results in the following formulation: Definition 3. Recent History Formulation for Causal Inference: for a pair of event labels (z, y), the binary treatment variable Zwt at time t is defined by whether or not z has occurred at least once within a window w into the past from t. The outcome variable is the occurrence rate of the effect label y at time t given the treatment, Yt = y|Zt(t). The covariates X w t at time t are a binary vector, depending on whether other event labels have occurred at least once in [t w, t), w 2 (0, t).
In event streams, it is possible for the cause and effect to be the same label in history, for instance, consider a system with negative self-feedback (eating too much sugar for example). Hence z and Xt could be identical to or contain y, but z and Xt must be disjoint sets of events.
We summarize the assumptions in this formulation of causal inference: 1) events before t w have no impact on y’s occurrence rate at time t. This enables memory in time yet provides a compact representation of history, 2) only the occurrence of z in the window impacts y’s rate at time t (regardless of the number of occurrences, hence we do not differentiate between multiple occurrences within w), and 3) the specific times of z’s occurrences do not further affect y’s rate at time t [16, 29, 15]. Such a model can be robust to outliers or noisy historical observations.

heading: 3.2 Defining Average Treatment Effect

text:
To measure how label y responds to historical occurrences of z, average treatment effect (ATE) [36] can be extended to our formulation in MPPs. We define the mean potential outcome Yt under treatment assignment Zwt = k as µky := 1 T R T 0 k y(t)dt, where T is the maximal time horizon. µky has an intuitive interpretation: it is the rate of y given Zwt = k in its history w, averaged over time horizon t0 = 0 to T . Note that µky measures the temporal average of instantaneous event occurrence counts under treatment. It is a natural way to aggregate an event’s expected occurrence counts in continuous time and captures the event dynamics affected by the proposed treatment. Definition 4. Average Treatment Effect (ATE) for event pairs is defined as:
ATE = EHT [µ 1 y µ 0 y] = EHT [ 1 T
Z
t
1y(t) 0 y(t)dt] (2)
where 1y(t) is y|Zt(t) at time t if z occurs at least once in Ht (Z w t = 1) and 0y(t) is y|Zt(t) at t if z does not occur in Ht (Zwt = 0). The EHT [·] is the expectation with respect to the random trajectories over [0, T ] as induced by the history dependent MPP. Note that since the rate ky is the instantaneous expected count of occurrences, ATE can also be viewed as a measure of expected count differences, normalized by the length of the time horizon T .
For a known multivariate point process, one can find the counterfactual ̂y|Zt(t) and hence the ATE by looking up the conditional intensity function with the corresponding parental conditions after making an adjustment to the treatment. Specifically, removing all occurrences of z in time period [t w, t) if the factual is Zwt = 1 or inserting z at time t w if the factual is Zwt = 0.
Theorem 1 provides justification for choosing ATE to measure pairwise causal relations in a MPP: Theorem 1. If z is not a direct cause of y in a multivariate point process, ATE for pair (z, y) is 0.
All detailed proofs are in the Appendices. The proof of Theorem 1 directly uses Definition 1 to show that the intensity rate does not change per z. The above theorem also implies that if z is not a parent of y in a graphical event model representation of the underlying MPP, the ATE for (z, y) is 0 as in Equation 2 from Definition 4.
To use ATE as defined in Equation 2, there are several assumptions that must hold in order to mimic the randomized trials to truly establish causal relationships. The ignorability condition indicates that whether y is 0 or 1 at each time t does not depend on whether Zwt = 1 or 0 at that time, i.e., ( 1y(t), 0y(t)) ? Zwt |Xt, 8t (although the y’s intensity rate still depends on Zwt ). Moreover, overlap condition states that each time t has a strictly positive chance for Zwt = 1 to happen given its history, i.e., 0 < P (Zwt = 1|Xt) < 1. Strong ignorability is often used when both overlap and ignorability hold true. If the assumption does not hold, possibly due to covariate differences in treatment groups, we need to adjust for these effects using propensity scores. This is particularly important in observational event data, where different times t with Zwt = 1 and Zwt = 0 may not be directly comparable, because the covariates Xwt may not be similar to each other. It can be studies further with the back-door criterion [31].

heading: 3.3 Defining Propensity Scores

text:
The propensity score is proposed to resolve covariate differences in non-randomized experiments to mimic a randomized study. It is a balancing score: conditioned on any balancing score, the distribution of observed covariates will be similar between the treated and control groups. It has been widely used in i.i.d. settings [22, 5, 13, 4].
We define a set of balancing scores, b(H zt ) (b⇤t for short), where the goal is to make treatment conditionally independent of covariates given the scores: Zt ? H Zt |b(H Z t ), 8t. Based on the specification provided by Definition 3, we use the recent window assumption and summarize the history H
z t with a set of covariates Xwt . Hence the balancing score would achieve Zwt ? Xwt |b(Xwt ), 8t. The most trivial b(Xwt ) is b(Xwt ) = Xwt . However, to consider the sampling distribution of potential cause Zwt , many-to-one functions of Xwt offer better adjustment, and the coarsest of such a function is the propensity score. We define X is finer than Y (equivalently, Y is coarser than X) if y = f(X). "Coarsest" here means that Y ’s dimension cannot be reduced further.
We use the above defined treatment Zwt and covariates Xwt and derive the equivalent form of propensity score for MPPs. The goal is to derive the propensity score e(Xwt ), or e⇤t for short, such that P (Zwt = 1|Xwt , e⇤t ) = P (Zwt = 1|e⇤t ). We know that:
P (Zwt = 1|X w t ) =
D(Zwt = 1;X w t )
D(Xwt ) , (3)
where D(Zwt = 1;Xwt ) is the duration that Zwt is observed true given observing covariates Xwt in the dataset. We define P (Zwt = 1|et) similarly as P (Zwt = 1|e⇤t ) = D(Zwt =1;e ⇤ t )
D(e⇤t ) . D(Zwt = 1; e⇤t )
can be interpreted as the duration that Zwt is observed true, i.e., the total duration for which z occurs in [t w, t], 8t, in the dataset with the same e⇤t in the relevant preceding windows. To ensure the required equality, we define e⇤t as: Definition 5. The propensity score for all covariates in history Xwt at time t is defined as e⇤t = e(Xwt ) = P (Z w t = 1|X w t ).
Note that Zwt at any time t follows a distribution P (Zwt = 1|Xwt ), which depends on the values of Xwt . This quantity can be estimated efficiently by counting the duration of Zwt = 1 with different Xwt values, which we will elaborate upon later. The major distinction between the proposed propensity score and that for the traditional i.i.d. case is that e⇤t is time-variant and incorporates all occurrences and non-occurrences of events in a window-viewed history.
We can show that times with the same value b⇤t but different treatment Zwt can act as controls for each other, as their expected outcome difference equals the proposed ATE: Theorem 2. Suppose the treatment event Zwt is strongly ignorable and b⇤t is a balancing score. Let [µky |A] := 1 T R t k y(t)|Adt. Then, [µ 1 y|X w t , Z w t = 1] [µ 0 y|X w t , Z w t = 0] = [µ 1 y|b ⇤ t ] [µ 0 y|b ⇤ t ].
Proof sketch. If the treatment is strongly ignorable, then it is also strongly ignorable given any balancing score. Therefore, we can show that the expected differences of both sides are equal. We further can show that treatment Zwt and covariates in history Xwt are conditionally independent, given the propensity score e⇤t , which makes e⇤t a balancing score. Please see the Appendices for exact proofs and further discussions.
It follows that the two-sampling process gives an unbiased estimate of ATE, using pair matching and sub-classification techniques [36] to adjust for the propensity score. However, in practice, pair
matching is difficult to perform given the continuous nature of time t which results in infinite sampling sizes. Subclassification is also difficult when the number of covariate event labels is large, leading to a high sample requirement on T .

heading: 4 Estimating ATE Scores

text:
There are several parameters that need to be provided or estimated in order to adjust for the propensity score and obtain ATE, including the window size for the treatment definition, conditional intensity rate of the treatment, and the outcome conditional or cumulative intensity rates.
Window Size w for the Treatment: The window size w of the treatment Zwt definition is provided as an input to the ATE estimation procedure. Moreover, it is used as the window size for all quantities associated with a window, namely Xwt , Zwt , y|Zt(t), and P (Z w t |X w t ). If one wants to estimate a window, practitioners should exercise prior knowledge to set w since the treatment effect with different windows will produce different results.
Alternately, one could take a data-driven approach to learning a window for a particular pair of event labels (z, y). One candidate is the heuristic pairwise approach in [7] that chooses from a set of inter-event times in the event stream data. Note that windows that are too short or too long could result in treatments that are always absent or present, respectively, making the ATE computation unsuitable and unstable. There is therefore a trade-off between long (which captures the historical influences as much as possible) and short windows (which captures only the important subset of history). We include more details in the Appendix.
Propensity Score e⇤t Estimation: Let D(Zwt = 1;xwt ) be the duration that Zwt is observed true in the dataset and that the condition xwt is true in the relevant preceding window w. Formally,
D(Zwt = 1;x w t ) =
N+1X
i=1
Z ti
ti 1
Ixwt (Z w t = 1, t)dt (4)
where Ixwt (Z w t = 1, t) is an indicator for whether each instantiation of Xwt , xwt , is true and Zwt = 1 at time t as a function of the relevant windows w. Similarly, D(xwt ) is the duration that the condition xwt is true in w. Hence, P (Zwt = 1|xwt ) = D(Zwt =1;x w t )
D(xwt ) and P (Zwt = 1|Xwt ) is a vector of such
probabilities with size |Xwt |.
Outcome Conditional Intensity Rate y|Zt(t): Given the most recent history view in the treatment definition, we estimate y|Zt(t), 80  t  T by using the same recent history formulation. This assumption has been used previously in MPPs, such as in proximal graphical event models (PGEMs) [7] where the rate of occurrence of an event type depends only on whether or not its parents have occurred in the most recent history.
When the parents U of all nodes V are known, the log likelihood of a MPP represented by a PGEM can be simplified to a function of counts and durations in the data and the conditional intensity rates:
logL(D) = X
V
X
u
v|uD(u) +N(v;u) ln( v|u)
(5)
where u is an instantiation of the parent of each variable v, N(v;u) is the number of times that X is observed in the dataset and that the condition u (from 2|U| possible parental combinations) is true in the relevant preceding windows, and D(u) is the duration over the entire time period where the condition u is true. Formally, N(v;u) = PN i=1 I(li = V )I wv u (ti) and D(u) = PN+1 i=1 R ti ti 1
Iwvu (t)dt, where Iwvu (t) is an indicator for whether u is true at time t as a function of the relevant windows wv . We refer the reader to [7] for more details around model learning.
We set z as the parent of y in a PGEM and then can estimate wy|z and w y|z̄ , which are intensities when z does and does not occur in the provided proximal window w. It is efficient since we only look at two events and only need to estimate intensity rates without learning windows. Once we have these two intensities, we can then look up y|Zt(t) by finding the treatment Z w t state at any query time t.
Outcome Cumulative Intensity Rate ⇤wfy|Zt(t): The cumulative intensity is computed as the integral of the conditional intensity over a future window wf , assuming that no other events occur. Again, given a learned point process model, we take a recent history view and compute the integral using
Algorithm 1 Inverse Probability Weighting for Events Input: Event data D, time horizon T Sample times t from 0 to T . Compute the propensity scores e⇤t , 8t, per Equation 3. Compute wt, 8t. per Equation 6 or stabilized scores. Train a PGEM My to estimate y|Zt(t), 8t. Predict with My for a given Ht the outcome 0y(t) with Zwt = 1 and 1y(t) with Zwt = 0. Compute ATE per Equation 7.
estimates for conditional intensities wy|z and w y|z̄ . Computing the integral from t to t + wf is straightforward as these are piece-wise constant conditional intensities. Furthermore, it is only possible for the treatment condition to change from Zwt = 1 to Zwt = 0 in the future window, when the condition from the occurrence of event label z changes, but not the other way around.
Inverse Probability Weighting for Events: The idea of weighting samples is as follows. If the time units for data involving the treatment Zwt = 1 group are far fewer than those in the control group Zwt = 0 in an event dataset, we need to upweight the time units in the treatment group and down-weight those in the control group for better estimation. We use inverse probability of treatment weighting (IPTW) as the weighting scheme, in order to adjust for this imbalance in the population.
We define the weight for the conditional intensity at time t in the treated Zwt = 1 group as ↵t = 1 P (Zwt =1|X w t ) = 1e⇤t , and for the control group Zwt = 0, the weight becomes ↵t = 1 1 P (Zwt =1|Ht)
= 1
1 e⇤t . Together, the weight for the outcome intensity rate at time t is defined as:
↵t = Zwt e⇤t + 1 Zwt 1 e⇤t w (6)
Using Equation 6, we propose an inverse probability weighting method for events, detailed in Algorithm 1. We first estimate the propensity score using Equation 3, and then estimate wt for all t using Equation 6. Then we can choose a PGEM My to predict y|Zt(t), for both the factual outcome and counterfactual y|Zt(t). Computing their empirical expectation provides an ATE estimate:
ATE = EHT [ 1 T
Z T
t=0
↵t · 1 y(t) 1 T
Z T
t=0
↵t · 0 y(t)] (7)
Equation 7 requires one to compare occurrence rates for all t in [0, T ], which includes some epochs with treatment and some without. Since this is a continuous time setting, integration over time is not straightforward. Hence, we propose using a sampling procedure to compute the inner integral over time in Equation 7 as 1S PS t=1 ↵t · 1 y(t) 1 S PS t=1 ↵t · 0 y(t), where S is the desired number of epoch samples from t0 = 0 to T . Corollary 3. Inverse probability weighting for events with balancing scores. Under strongly ignorable treatment, the weighted average with inverse probability weighting is unbiased for ATE when
weights equal those in Equation 6.
Proof Sketch. The inner integral of ATE can be shown to equal to [µ1y] [µ0y], and further [µ1y(t)] = 1 T R t ZWt · ⇤ y(t) e⇤t dt and [µ0y(t)] = 1 T R t (1 ZWt )· ⇤ y(t) 1 e⇤t dt.
We show in Corollary 3 that our ATE definition is an unbiased estimator. Consistency results, however, rely on the estimation results of intensities from specific point process models. For example, when the window sizes are given, PGEM estimator [7] for intensities are consistent.
Since we consider (interval) stationary point process here, we can follow the standard computational procedure to estimate the sample variance of ATE estimator and a sample-based confidence interval for the ATE estimate (for example, results from Theorem 6.2 in [19] can be directly extended to our setting), as the samples are from the same distribution.
Factual and Counterfactual Treatment Usages To compute the ATE, we need both a factual and counterfactual treatment for time. For time t with factual Zwt = 1, the counterfactual Zwt = 0 treatment is modeled by removing all occurrences of z in [t w, t). For time t with factual Zwt = 0, the counterfactual Zwt = 1 treatment is modeled by introducing one occurrence of z in [t w, t).
(a) ATE accuracy results on different causal pair scores for 3 PGEM synthetic datasets.
(b) ATE accuracy results on different causal pair scores for 3 Hawkes synthetic datasets.
For ease of computation, we introduce the counterfactual occurrence of z at t w for instantaneous conditional intensity computation, although its exact location would not impact the outcome due to the proximal assumption. For cumulative intensity outcomes, we introduce the counterfactual occurrences at t t, with t = w2 which is the middle point of the window. Other choices t of introducing the counterfactual treatment are also possible.
Stabilized IPTW One common issue with IPTW is that sometimes the propensity scores for some time unit t can get very close to 0, which indicates Z is extremely unlikely to occur in the window [t w, t); this is even more likely for continuous time data. Hence, the weights for those t become extremely large, causing unstable estimations. To combat this issue, the stabilized IPTW [35] uses the marginal probability of treatment to counteract such an instability. It is formulated as ↵t = Zwt ·P (Z w t =1)
e⇤t + (1 Z
w t )·(1 P (Z w t =1))
1 e⇤t .

heading: 5 Empirical Evaluation

text:
Evaluating causal inference algorithms is more difficult than those for prediction tasks since observational datasets rarely contain the ground truth treatment effects. To this end, most experiments in the literature analyze causal models using a synthetic dataset where the ground truth is known [20]. We begin by comparing the ATE estimation performance of our proposed IPTW methods on synthetic event datasets, generated using different parameters. Per standard practice in causal inference literature, we use root mean squared error (RMSE) to measure the ATE accuracy of each method, along with its standard deviation. We first focus on the experiments with the instantaneous conditional intensity outcome. All experiments are done on a machine with 2.9 GHz quad-core CPU.
Since outcome y is not directly observable, simple adaptation of ATE from the i.i.d. case would not work. There are not many well-established baselines that use event occurrence rate as the outcome in multivariate point processes for comparison. We adapt two baseline scores from fitting parametric models to the intensity rates from cause-effect associate scores 11[8]: CI (conditional intensity) and CIM , which consider a single parent event and a set, respectively. For the first baseline, we consider the association between a pair of events, (z, y) and assume that the intensity of y only depends on whether or not z has occurred at least once within a specified time window w. Hence we define a conditional intensity score to estimate the causal effect of z on y as CI(z, y) = wy|z w y|z̄ . We extend this definition to allow for y depending on the historical arrival of a set of events instead of
just a single z. This means that there are possibly 2|X|+1 conditional intensity rates, and we need to aggregate the score for z given all other conditions of the parent set. For our purposes we will use the mean of y over all settings. Formally, CIM (z, y) = mean ⇣ wy|z,x w y|z̄,x ⌘ , when z is a parent of y and otherwise the score is 0.
We compare CI scores with three versions of ATE estimation based on our proposed IPTW methods. First, we compute ATE with no weighting as per Equation 2; labeled as IP-NW. Second, we use the proposed IPTW with non-stable weight (IP-NS) with weights as per Equation 6. Last, we use IPTW with stable weight (IP-Stable) to compute ATE.
Synthetic PGEM Datasets To test our estimation methods on synthetic data, we first generate event data that adheres to the proximal assumption for the intensity functions [7]. We generate 3 models with different numbers of events, randomly generated graph structures among events, fixed window size of w = 30, T = 2000, and random intensities between 0.1 and 0.4. We use the data and the generated model to obtain the true estimates of y|Zt(t) at chosen times t and hence can compute the ground truth ATE. Since we observed that the sample size S of t (103 to 105) in the ATE estimation does not impact the results much, we use sample size S = 103 for all our experiments. Shown in
Figure 1a on three synthetic PGEM datasets, non-weighted IP-NW performs worse than CI scores when window size w is small but better with large windows. IP-Stable’s RMSE is consistently the smallest across datasets and window sizes, sometimes with a 6 times reduction in error from CI scores. We did not include IP-NS in this figure as it has large errors in comparison to other methods. IP-NS with non-stable weight can lead to non-stable estimation. For exact numerical values, please refer to the appendix.
Synthetic Hawkes Dataset We also test our approaches on synthetic multivariate Hawkes process datasets using an existing toolbox [6]. We again generate 3 datasets with 30, 40, and 50 event labels, and w = 15. We use a fixed base rate 0.016 and each parental event leads to additive spike of 0.06 to the base rate, with an exponential decay rate 0.15. We generate event streams with T = 2000. For estimating the counterfactual rate, ̂y|Zt(t), to compute ATE in the case where there is no treatment event over [t w, t), we introduce a counterfactual treatment event at time t w. Figure 1b shows the RMSE results of different algorithms. IP-NW outperforms CI scores in all cases. IP-Stable shows the lowest RMSE among all methods, generally achieving 3 to 10 times better accuracy. In addition, the performance of IP-NW and IP-Stable are relatively stable with respect to the window sizes, with minimal changes. For exact numerical values, please refer to the appendix.
Hybrid Dataset We also generate a synthetic hybrid dataset that combines proximal graphical event models with additive excitation similar to Hawkes processes with a constant kernel. For further details about the generating process and the results, please refer to Appendices G and H. IP-NW outperforms CI scores in all cases but one, and IP-Stable shows the lowest RMSE in all but two cases.
Diabetes Dataset We also test our methods on the diabetes dataset [14] – a real-world dataset which we process into events for meals, exercise activity, insulin dosage and changes in blood glucose measurements for 70 diabetes patients. We treat the assessments in [2] as the ground truth, where an expert provided 11 pairs such that a cause label is more likely to make the effect label occur. Since the assessments are only partial and do not provide the true ATE, we use hits@K among highest absolute estimated ATE values to measure performance in this experiment, which is a popular metric for information retrieval. Specifically, we determine how many of the 11 pairs are recovered by a method’s top K absolute scores. The dataset is split into 50%/50% training/test sets, and optimal window setting is determined on the training set, which is then deployed in the test set for evaluation. w = {0.1, 0.3, 0.5, 1} days for all models were considered during training. Table 1 compares the Hits@K across all methods on the test set for K = {10, 15, 20}, illustrating that IP-Stable is able to recover more of the expert-assessed pairs in two out of three K values.
Cumulative Intensity Outcome Experiments We also test our framework with the outcome being the cumulative intensity. The procedure can be done directly by replacing y(t) with its cumulative version ⇤wfy (t), computed from a future window [t, t+ wf ), as discussed in Section 4. Hence the overall computation becomes ATE = EHT [ 1 T R T t=0 ↵t · ⇤ 1 y(t) 1 T R T t=0 ↵t · ⇤ 0 y(t)] . To adapt the CI baselines to cumulative setting, we compute the cumulative version of CI, CCI , with CCI(t) = CIt ⇥ wf . We also name our methods CIP as the cumulative version of the IP methods. All other computations remain the same. We are given the future window size wf = 15 for all methods. To explore the data-driven approach for learning the window, we do not provide w here and instead let it learn from data [7]. We show the results of ATE for the cumulative intensity outcome with 2 PGEMs in Table 2. The results again confirm better performance of CIP methods over baselines, specially NW and Stable versions. CIP-NS can be unstable as indicated previously.

heading: 6 Conclusion

text:
We have proposed a framework for pairwise event causality in a multivariate point process and formalized the problem with novel definitions of treatment, outcome, and propensity scores. Our definitions allow for efficient modeling of data where outcomes occur multiple times on the timeline. We estimated the average treatment effect using a propensity score weighting procedure that achieves the best performance against baselines. Our work bridges causal inference with multivariate point processes, showing promising performance in estimating pairwise causal relationships among events. Our framework can be extended in many ways. For example, ATE with multiple causes could be computed if we specify an outcome model given such multiple causes. Our treatment and covariate definitions focus on binary values but could be extended to marked point processes where the event has a real valued measurement, by choosing an appropriate aggregation function and outcome models. Future work could study efficient estimation approaches of ATE without sampling and different historical representations of the treatment. Another direction is a count-based estimand over the current proximal assumption, which would result in multiple treatment value problems.
Limitations And Societal Impact: Our method focuses on potentially long event sequences with multiple occurrences of all events. We assume proximal historical influence, along with strong ignorability and no confounding. The usage of causal models should be cautioned, to draw potentially harmful conclusions. In particular, the window sizes in our formulation, which can be subject to human inputs, can induce biases and provide different conclusions with different values.

heading: Acknowledgement

text: We thank to anonymous reviewers who have provide helpful comments. This research was partly developed with funding from the Defense Advanced Research Projects Agency (DARPA), under Program No. FA8750-19-C-0206. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies, either expressed or implied, of DARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for governmental purposes notwithstanding any copyright annotation therein 1. XS is supported by supported by IBM AI Horizons Network. Nicholas Mattei was supported by an IBM Academic Award.

title: Causal Inference for Event Pairs in Multivariate Point Processes


summaries: This paper discusses how message passing and feature aggregation are affected by abnormal node features, and designs a message passing module called Adaptive Message Passing (AMP) with adaptive residual connection and feature aggregation. In this way, the module can adaptively deal with the node that is inconsistent with the local feature during the message passing. The paper proposes the AirGNN that is built up with AMP as the message-passing layer in GNN, the experiment result shows that the AirGNN gets a more robust result on both the abnormal node and normal node classification comparing with the baselines.

summaries: This paper targets to design GNNs with stronger resilience to abnormal node features. Empirical examination is firstly conducted by replacing the features of randomly selected nodes with random Gaussian noise. The comparison between models with and without the residual connections show the helpfulness of residual connections on learning with abnormal features.

summaries: This paper observed interesting phenomenons: feature aggregation helps smooth out abnormal features but would lead to over-smoothing for normal features, residual link helps adjust feature smoothness for normal features but may hurt the performance when we have abnormal features. Therefore, it points out that we need to tradeoff between normal and abnormal features while designing GNNs. Then, based on the observation, this work proposed the AMP method and the AirGNN model, which are shown by experiments to be effective under various abnormal feature scenarios.

summaries: This paper examines the behavior of GNNs when node features are perturbed during test time.  AirGNN is proposed based on the observations from the standard GNNs. AirGNN adaptively computes the weight of the initial embedding and aggregated embeddings and is robust to random/adversarial noise keeping the standard classification performance.

abstractText: Graph neural networks (GNNs) have shown the power in graph representation learning for numerous tasks. In this work, we discover an interesting phenomenon that although residual connections in the message passing of GNNs help improve the performance, they immensely amplify GNNs’ vulnerability against abnormal node features. This is undesirable because in real-world applications, node features in graphs could often be abnormal such as being naturally noisy or adversarially manipulated. We analyze possible reasons to understand this phenomenon and aim to design GNNs with stronger resilience to abnormal features. Our understandings motivate us to propose and derive a simple, efficient, interpretable, and adaptive message passing scheme, leading to a novel GNN with Adaptive residual, AirGNN1. Extensive experiments under various abnormal feature scenarios demonstrate the effectiveness of the proposed algorithm.

authors:
Yao Ma
Jiliang Tang

title: Deep Learning on Graphs

authors:
Justin Gilmer
Samuel S Schoenholz
Patrick F Riley
Oriol Vinyals
George E Dahl

title: Neural message passing for quantum chemistry

venue: In International Conference on Machine Learning,

authors:
Thomas N Kipf
Max Welling

title: Semi-supervised classification with graph convolutional networks

venue: arXiv preprint arXiv:1609.02907,

authors:
Michaël Defferrard
Xavier Bresson
Pierre Vandergheynst

title: Convolutional neural networks on graphs with fast localized spectral filtering

venue: In Proceedings of the 30th International Conference on Neural Information Processing Systems,

authors:
William L Hamilton
Rex Ying
Jure Leskovec

title: Inductive representation learning on large graphs

venue: arXiv preprint arXiv:1706.02216,

authors:
Franco Scarselli
Marco Gori
Ah Chung Tsoi
Markus Hagenbuchner
Gabriele Monfardini

title: The graph neural network model

venue: IEEE transactions on neural networks,

authors:
Zonghan Wu
Shirui Pan
Fengwen Chen
Guodong Long
Chengqi Zhang
S Yu Philip

title: A comprehensive survey on graph neural networks

venue: IEEE transactions on neural networks and learning systems,

authors:
Ming Chen
Zhewei Wei
Zengfeng Huang
Bolin Ding
Yaliang Li

title: Simple and deep graph convolutional networks

venue: Proceedings of the 37th International Conference on Machine Learning,

authors:
Johannes Klicpera
Aleksandar Bojchevski
Stephan Günnemann

title: Predict then propagate: Graph neural networks meet personalized pagerank

venue: arXiv preprint arXiv:1810.05997,

authors:
Yao Ma
Xiaorui Liu
Tong Zhao
Yozen Liu
Jiliang Tang
Neil Shah

title: A unified view on graph neural networks as graph signal denoising

venue: Proceedings of the 30th ACM International Conference on Information and Knowledge Management,

authors:
Hoang Nt
Takanori Maehara

title: Revisiting graph neural networks: All we have is low-pass filters

venue: arXiv preprint arXiv:1905.09550,

authors:
Lingxiao Zhao
Leman Akoglu

title: Pairnorm: Tackling oversmoothing in gnns

venue: In International Conference on Learning Representations,

authors:
Qimai Li
Zhichao Han
Xiao-Ming Wu

title: Deeper insights into graph convolutional networks for semi-supervised learning

venue: In Proceedings of the AAAI Conference on Artificial Intelligence,

authors:
Kenta Oono
Taiji Suzuki

title: Graph neural networks exponentially lose expressive power for node classification

authors:
Miller McPherson
Lynn Smith-Lovin
James M Cook

title: Birds of a feather: Homophily in social networks

venue: Annual review of sociology,

authors:
Dengyong Zhou
Olivier Bousquet
Thomas Navin Lal
Jason Weston
Bernhard Schölkopf

title: Learning with local and global consistency

authors:
Meng Liu
Hongyang Gao
Shuiwang Ji

title: Towards deeper graph neural networks

venue: In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. ACM,

authors:
Prithviraj Sen
Galileo Namata
Mustafa Bilgic
Lise Getoor
Brian Galligher
Tina Eliassi- Rad

title: Collective classification in network data

venue: AI magazine,

authors:
Oleksandr Shchur
Maximilian Mumme
Aleksandar Bojchevski
Stephan Günnemann

title: Pitfalls of graph neural network evaluation

venue: arXiv preprint arXiv:1811.05868,

title: Open graph benchmark: Datasets for machine learning on graphs

authors:
Dingyuan Zhu
Ziwei Zhang
Peng Cui
Wenwu Zhu

title: Robust graph convolutional networks against adversarial attacks

venue: In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining,

authors:
Diederik P Kingma
Jimmy Ba

title: Adam: A method for stochastic optimization

venue: arXiv preprint arXiv:1412.6980,

authors:
Daniel Zügner
Amir Akbarnejad
Stephan Günnemann

title: Adversarial attacks on neural networks for graph data

venue: In KDD

authors:
Yaxin Li
Wei Jin
Han Xu
Jiliang Tang

title: Deeprobust: A pytorch library for adversarial attacks and defenses, 2020

authors:
Xuran Pan
Song Shiji
Huang Gao

title: A unified framework for convolution-based graph neural networks

authors:
Meiqi Zhu
Xiao Wang
Chuan Shi
Houye Ji
Peng Cui

title: Interpreting and unifying graph neural networks with an optimization framework, 2021

authors:
Siheng Chen
Yonina C. Eldar
Lingxiao Zhao

title: Graph unrolling networks: Interpretable neural networks for graph signal denoising, 2020

authors:
Guohao Li
Matthias Müller
Ali Thabet
Bernard Ghanem

title: Deepgcns: Can gcns go as deep as cnns

venue: In The IEEE International Conference on Computer Vision (ICCV),

authors:
Keyulu Xu
Chengtao Li
Yonglong Tian
Tomohiro Sonobe
Ken-ichi Kawarabayashi
Stefanie Jegelka

title: Representation learning on graphs with jumping knowledge networks

venue: Proceedings of the 35th International Conference on Machine Learning,

authors:
Negin Entezari
Saba A Al-Sayouri
Amirali Darvishzadeh
Evangelos E Papalexakis

title: All you need is low (rank) defending against adversarial attacks on graphs

venue: In Proceedings of the 13th International Conference on Web Search and Data Mining,

authors:
Wei Jin
Yao Ma
Xiaorui Liu
Xianfeng Tang
Suhang Wang
Jiliang Tang

title: Graph structure learning for robust graph neural networks

venue: In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining,

authors:
Yu Chen
Lingfei Wu
Mohammed Zaki

title: Iterative deep graph learning for graph neural networks: Better and robust node embeddings

venue: Advances in Neural Information Processing Systems,

authors:
Xiaorui Liu
Wei Jin
Yao Ma
Yaxin Li
Hua Liu
Yiqi Wang
Ming Yan
Jiliang Tang

title: Elastic graph neural networks

venue: Proceedings of the 38th International Conference on Machine Learning,

authors:
Wei Jin
Yaxing Li
Han Xu
Yiqi Wang
Shuiwang Ji
Charu Aggarwal
Jiliang Tang

title: Adversarial attacks and defenses on graphs

venue: ACM SIGKDD Explorations Newsletter,

authors:
Yanqiao Zhu
Weizhi Xu
Jinghao Zhang
Qiang Liu
Shu Wu
Liang Wang

title: Deep graph structure learning for robust representations: A survey

venue: arXiv preprint arXiv:2103.03036,

title: A primal–dual splitting method for convex optimization involving lipschitzian, proximable and linear composite terms

venue: Journal of optimization theory and applications,

heading: 1 Introduction

text:
Recent years have witnessed the great success of graph neural networks (GNNs) in representation learning for graph structure data [1]. Essentially, GNNs generalize deep neural networks (DNNs) from regular grids, such as image, video and text, to irregular data such as social, energy, transportation, citation, and biological networks. Such data can be naturally represented as graphs with nodes and edges. The key building block for such generalization is the neural message passing framework [2]:
x(k+1)u = UPDATE (k) ( x(k)u ,m (k) N (u) ) (1)
where x(k)u ∈ Rd denotes the feature vector of node u in the k-th iteration of message passing, and m
(k) N (u) is the message aggregated from u’s neighborhood N (u). The specific design of message passing scheme can be motivated from spectral domain [3, 4] or spatial domain [5, 6, 7, 2]. It usually linearly smooths the features in a local neighborhood on the graph.
GNNs have achieved superior performance in a large number of benchmark datasets [8] where the node features are assumed to be complete and informative. However, in real-world applications, some node features could be abnormal from various aspects. For instance, in social networks, new users might not have complete profile before they make connections with others, leading to missing user features. In transportation networks, node features can be noisy since there exist certain
1The implementation is available at https://github.com/lxiaorui/AirGNN.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
uncertainty and dynamics in the observation of the traffic information. What is worse, node features can be adversarially chosen by the attacker to maliciously manipulate the prediction made by GNNs. Therefore, it is greatly desired to design GNN models with stronger resilience to abnormal node features.
In this work, we first perform empirical investigations on how representative GNN models behave on graphs with abnormal features. Specifically, based upon standard benchmark datasets, we simulate the abnormal features by replacing the features of randomly selected nodes with random Gaussian noise. Then the performance of node classification on abnormal features and normal features are examined separately. From our preliminary study in Section 2, we reveal two interesting observations: (1) Feature aggregation can boost the resilience to abnormal features, but too many aggregations could hurt the performance on both normal and abnormal features; and (2) Residual connection helps GNNs benefit from more layers for normal features, while making GNNs more fragile to abnormal features. We then provide possible explanations to understand these observed phenomena from the perspective of graph Laplacian smoothing. Our analyses imply that there might exist an intrinsic tension between feature aggregation and residual connection, which results in a performance tradeoff between normal features and abnormal features.
Motivated by these findings and understandings, we aim to design new GNNs with stronger resilience to abnormal features while largely maintaining the performance on normal features. Our contributions can be summarized as follows:
• We discover an intrinsic tension between feature aggregation and residual connection in GNNs, and the corresponding performance tradeoff between abnormal and normal features. We also analyze possible reasons to explain and understand these findings.
• We propose a simple, efficient, principled and adaptive message passing scheme, which leads to a novel GNN model with adaptive residual, named as AirGNN.
• Extensive experiments under various abnormal feature scenarios demonstrate the superiority of the proposed algorithm. The ablation study demonstrates how the adaptive residuals mitigate the impact of abnormal features.

heading: 2 Preliminary

text:
Before introducing the preliminary study, we first define the notations used throughout the paper.
Notations. We use bold upper-case letters such as X to denote matrices. Given a matrix X ∈ Rn×d, we use Xi to denote its i-th row and Xij to denote its element in i-th row and j-th column. The
Frobenius norm and `21 norm of a matrix X are defined as ‖X‖F = √∑
ijX 2 ij and ‖X‖21 =∑ i ‖Xi‖2 = ∑ i √∑ jX 2 ij , respectively. We define ‖X‖2 = σmax(X) where σmax(X) is the
largest singular value of X.
Let G = {V, E} be a graph with the node set V = {v1, . . . , vn} and the undirected edge set E = {e1, . . . , em}. We use N (vi) to denote the neighboring nodes of node vi, including vi itself. Suppose that each node is associated with a d-dimensional feature vector, and the features for all nodes are denoted as Xfea ∈ Rn×d. The graph structure G can be represented as an adjacent matrix A ∈ Rn×n, where Aij = 1 when there exists an edge between nodes vi and vj , and Aij = 0 otherwise. The graph Laplacian matrix is defined as L = D − A, where D is the diagonal degree matrix. Let us denote the commonly used feature aggregation matrix in GNNs [3] as Ã = D̂− 1 2 ÂD̂− 1 2 where Â = A+ I is the adjacent matrix with self-loop and its degree matrix is
D̂. The corresponding Laplacian matrix is defined as L̃ = I− Ã. In this work, we focus on the setting where a subset of nodes in the graph contain abnormal features, while the remaining nodes have normal features. In the remaining of this paper, we use abnormal/normal features to denote nodes with abnormal/normal features, for simplicity.

heading: 2.1 Preliminary Study

text:
Experimental setup. To investigate how GNNs behave on abnormal and normal node features, we design semi-supervised node classification experiments on three common datasets (i.e., Cora,
CiteSeer and PubMed), following the data splits in the work [3]. Moreover, we simulate the abnormal features by assigning 10% of the nodes with random features sampled from a standard Gaussian distribution. The experiments are performed on representative GNN models covering coupled and decoupled architectures, including GCN [3], GCNII [9], APPNP [10], and their variants with or without residual connections in feature aggregations, denoted as w/Res and wo/Res. All methods follow the hyperparameter settings in their original papers. We examine how these models perform when the number of layers increases. Note that for the decoupled architectures such as APPNP, we fix the 2-layer MLP and increase the number of propagation layers. While for the coupled architectures such as GCN and GCNII, we increase the number of feature transformation and propagation layers simultaneously. We report the average performance over 10 times of random selection of the noise node sets. The node classification accuracy (mean and standard variance) on nodes with abnormal and normal features is illustrated in Figure 1 and Figure. 2, separately.
Observations. From Figure 1 and Figure 2, we can make the following observations: (1) Without residual connection, more layers (e.g., > 2 for GCN and GCNII, > 10 for APPNP) hurt the accuracy on nodes with normal features. However, more layers boost the accuracy on nodes with abnormal features significantly, before finally starting to decrease; (2) With residual connection, the accuracy on nodes with normal features keeps increasing with more layers2. However, the accuracy on nodes with abnormal features only increases marginally when stacking more layers, and then starts to decrease. While we only present the experiments on Cora, we defer the results on other datasets to Appendix C, which provide similar observations. To conclude, we can summarize these observations into two major findings:
• Finding I: Feature aggregation can boost the resilience to abnormal features, but too many aggregations could hurt the performance on both normal and abnormal nodes;
• Finding II: Residual connection helps GNNs benefit from more layers for nodes with normal features, while making GNNs more fragile to abnormal features.
2GCN w/Res is an exception because its residual is not appropriate, which is consistent with the experiments in the work [3].

heading: 2.2 Understandings

text:
In this subsection, we provide the understanding and explanation for aforementioned findings, from the perspective of graph Laplacian smoothing.
Understanding Finding I: Feature aggregation as Laplacian smoothing
The message passing in GCN [3], GCNII wo/ residual and APPNP wo/ residual (as well as many popular GNN models), follows the feature aggregation
Xout = ÃXin, (2)
where Xin and Xout represent the features before and after message passing layer, respectively. It can be interpreted as one gradient descent step for the Laplacian smoothing problem [11]
argmin X∈Rn×d
L1(X) := 1 2 tr ( X>(I− Ã)X ) = 1 2 ∑ (vi,vj)∈E ‖ Xi√ di + 1 − Xj√ dj + 1 ‖22, (3)
where di is the node degree of node vi. Eq. (2) can be derived from Xout = Xin−(I−Ã)Xin = ÃXin, with the initialization X = Xin and stepsize γ = 1. The Laplacian smoothing problem penalizes the feature difference between neighboring nodes. To reduce this penalty, the feature aggregation in Eq. (2) smooths the node features by taking the average of local neighbors, and thus can be considered as low-pass filter which gradually filters out high-frequency signals [12, 13]. Therefore, it increases the resilience to abnormal features which are likely to be high-frequency signals. In other words, the local neighboring nodes help to correct the abnormal features. Unfortunately, if applied too many times, these low-pass filters could overly smooth the features (well-known as oversmoothing [14, 15]) such that nodes are not distinguishable enough, providing an explanation to the degraded performance on both abnormal and normal features when stacking too many layers.
Understanding Finding II: Residual connection maintains feature proximity
To adjust the feature smoothness for better performance, APPNP [10] utilizes residual connections in message passing as follows
Xk+1 = (1− α)ÃXk + αXin, (4)
where X0 = Xin. It can be considered as an iterative solution for the regularized Laplacian smoothing problem [11]
argmin X∈Rn×d
L2(X) := α
2(1− α) ‖X−Xin‖2F +
1 2 tr ( X>(I− Ã)X ) , (5)
with initialization X = Xin and stepsize γ = 1− α due to Xk+1 = Xk − (1− α) ( α 1− α (Xk −Xin) + (I− Ã)Xk ) = (1− α)ÃXk + αXin.
GCNII [9] adopts a similar message passing but further combines a feature transformation layer in each message passing step, which leads to a coupled architecture, as contrast to the decoupled architecture of APPNP. The residual connection naturally arises when regularizing the proximity between input and output features, as showed in the first term of L2(X). Such proximity can help avoid the trivial solution for the problem in Eq. (3), i.e., totally oversmoothed features only depending on node degrees, and consequently mitigates the oversmoothing issue. More intuitively, residual connections in GNNs provide direct information flows between layers that can preserve some necessary high-frequency signals for better discrimination between classes. More layers with residual provide a more accurate solution to Eq. (5), which explains the performance gain from deeper GNNs. Unfortunately, these residual connections also undesirably carry on abnormal features which are detrimental, leading to the inferior performance on abnormal features.

heading: 3 The Proposed Framework

text: In this section, we first motivate the proposed adaptive message passing scheme (AMP) with further discussions on our preliminary study. We then introduce more details about AMP, its interpretations, convergence guarantee and computation complexity, as well as the model architecture of AirGNN.

heading: 3.1 Design Motivation

text:
Our preliminary study in Section 2 reveals an intrinsic tension between feature aggregation and residual connection: (1) feature aggregation helps smooth out abnormal features, while it could cause inappropriate smoothing for normal features; (2) residual connection is essential for adjusting the feature smoothness, but it could be detrimental for abnormal features. Although this conflict can be partially mitigated by adjusting the residual connection such as the residual weight α in GCNII [9] and APPNP [10], such global adjustment cannot be adaptive to a subset of the nodes, e.g., the nodes with abnormal features. This is crucial because in practice we often encounter the scenario where only a subset of nodes contain abnormal features. Therefore, how to reconcile this dilemma still desires dedicated efforts. We then naturally ask a question: Can we design a better message passing scheme with node-wise adaptive feature aggregation and residual connection?
The motivation of the proposed idea builds upon the following intuition: while it is important to maintain the proximity between input and output features as in Eq. (5), it could be over aggressive to penalize their deviations by the square of Frobenius norm, i.e., ‖X−Xin‖2F = ∑n i=1 ‖Xi− (Xin)i‖22. The fact that this penalty does not tolerate large deviations weakens the capability to remove abnormal features through Laplacian smoothing. This motivates us to consider an alternative proximity penalty
‖X−Xin‖21 := n∑ i=1 ‖Xi − (Xin)i‖2, (6)
which instead penalizes the deviations by the `1 norm of row-wise `2 norms, namely `21 norm. The `21 norm promotes row sparsity in X−Xin, and it also allows large deviations because the penalty on large values is less aggressive, leading to the potential removal of abnormal features. Therefore, we propose the following Laplacian smoothing problem regularized by `21 norm proximity control:
argmin X∈Rn×d
λ‖X−Xin‖21 + 1
2 tr(X>(I− Ã)X), (7)
where λ ∈ [0,∞) is a parameter to adjust the balance between proximity and Laplacian smoothing. In order to easy the tuning of λ, we made a modification of Eq. (7):
argmin X∈Rn×d
L(X) := λ‖X−Xin‖21 + (1− λ)tr(X>(I− Ã)X), (8)
where λ ∈ [0, 1] controls the balance.

heading: 3.2 Adaptive Message Passing

text:
L(X) is a composite objective with non-smooth and smooth components. We optimize it by proximal gradient descent [16] and obtain the following iterations as the adaptive message passing (AMP):
Yk = Xk − 2γ(1− λ)(I− Ã)Xk = ( 1− 2γ(1− λ) ) Xk + 2γ(1− λ)ÃXk (9)
Xk+1 = argmin X
{ λ‖X−Xin‖21 + 1
2γ ‖X−Yk‖2F
} (10)
where X0 = Xin and γ is the stepsize to be specified later. Let Z = X−Xin, and Eq. (10) can be rewritten as:
Zk+1 = argmin Z
{ λ‖Z‖21 + 1
2γ ‖Z− (Yk −Xin)‖2F } = proxγλ‖·‖21(Y k −Xin) (11)
Xk+1 = Xin + Z k+1. (12)
The i-th row of the proximal operator in Eq. (11) can be computed analytically( proxγλ‖·‖21(X) ) i =
Xi ‖Xi‖2 max(‖Xi‖2 − γλ, 0) = max(1− γλ ‖Xi‖2 , 0) ·Xi. (13)
Note that the proximal operator returns 0 if the input vector is 0. Substituting X in Eq. (13) with Yk −Xin and combining Eq. (11) and Eq. (12), then Eq. (12) becomes
Xk+1i = (Xin)i + βi(Y k i − (Xin)i) = (1− βi)(Xin)i + βiYki , ∀i ∈ [n], (14)
where βi := max(1 − γλ‖Yki −(Xin)i‖2 , 0). To summarize, the proposed adaptive message passing (AMP) scheme is showed in Figure 4, and a diagram is showed in Figure 3. In detail, AMP works as follows:
• The first step takes a feature aggregation within the local neighbors with a self-loop weighted by 1− 2γ(1− λ);
• The second step computes a weight βi ∈ [0, 1] for each node vi depending on the local deviation ‖Yki − (Xin)i‖2.
• The final step takes a linear combination of input features Xin and the aggregated features Yk, where the node-wise residual is adaptively weighted by 1− βi for each node vi.
The convergence guarantee of AMP and parameter setting for the stepsize γ are illustrated in Theorem 1 and proved in Appendix A. According to Theorem 1, if we set γ = 14(1−λ) or γ = 1 2(1−λ) , then the first step of AMP can be simplified as Yk = 12X k + 12ÃX
k and Yk = ÃXk, respectively. The choice of stepsize will only impact the convergence speed but not the ultimate effect of AMP when it convergences to the fixed point solution. We also discuss the computation complexity per iteration of AMP in Remark 1. Theorem 1 (Convergence of AMP). Under the stepsize setting γ < 1
(1−λ)‖L̃‖2 , the proposed adaptive
message passing scheme (AMP) in Eq. (9) and Eq. (10) converges to the optimal solution of the problem defined in Eq. (8). In practice, it is sufficient to choose any γ < 12(1−λ) since ‖L̃‖2 ≤ 2. Moreover, if the connected components of the graph G are not bipartite graphs, it is sufficient to choose γ = 12(1−λ) since ‖L̃‖2 < 2. Remark 1 (Computation complexity). AMP is as efficient as simple feature aggregation Xout = ÃXin because the additional computation cost from the second and third steps in Figure 4 is in the order O(nd), where n is the number of nodes and d is the feature dimension. This is negligible compared with the computation cost O(md) in feature aggregation, where m is the number of edges, due to the fact that usually there are many more edges than nodes in real-world graphs, i.e., m n.

heading: 3.3 Interpretation of AMP

text:
Interestingly, the proposed AMP has a simple and intuitive interpretation as adaptive residual connection, which aligns well with our design motivation:
• If the feature of node vi, i.e., (Xin)i, is significantly inconsistent with its local neighbors, i.e., the aggregated feature Yki , then the local deviation ‖Yki − (Xin)i‖2 will be large, which leads to a βi close to 1. Therefore, the final step will assign a small weight to the residual, i.e., (1− βi)(Xin)i, and the aggregated feature Yki will dominate.
• On the contrary, if (Xin)i is already consistent with its local neighbors, ‖Yki − (Xin)i‖2 will be small, which leads to a βi close to 0. Thus, the residual will dominate, which is reasonable since there is less need to aggregate features in this case.
• To summarize, the local deviation ‖Yki − (Xin)i‖2 provides a natural transition from βi → 1 to βi → 0, and the transition can be modulated by λ which can be either learned or tuned as a hyperparameter through cross-validation. This transition provides an node-wise adaptive residual connection for the message passing scheme.
Adaptivity for abnormal & normal features. According to the homophily assumption on graph structure data [17, 18, 19, 3], the feature representations of normal features should be more consistent with local neighbors than abnormal features. As a result, AMP will assign more residual (i.e., smaller β) to normal features but less residual (i.e., larger β) to abnormal features, providing a customized tradeoff between feature aggregation and residual connection. Consequently, it can promote both the resilience to abnormal features and the performance on normal features. Above discussion also implies a clear physical meaning for β in AMP, and we formally define it as the adaptive score. Definition 1 (Adaptive score). The variables {β1, · · · , βn} in the adaptive message passing scheme (AMP) are defined as the adaptive scores for nodes {v1, · · · , vn} respectively in graph G. In particular, the larger βi is, the more likely the feature of node vi is abnormal. Remark 2 (Nonlinear smoother). Different from most existing message passing scheme which are linear smoothers, AMP is a nonlinear smoother because the weights {βi} are computed from Yk and Xin. This nonlinearity is the key to achieve adaptive residual connection for different nodes.

heading: 3.4 The Model Architecture

text:
The proposed adaptive message passing (AMP) can be used as a building block in many GNN models to improve the resilience to abnormal node features. In this work, we choose the the decoupled architectures as APPNP [10] and DAGNN [20], and propose the Adaptive residual GNN (AirGNN):
Xin = hθ(Xfea), (15) Ypre = AMP (Xin,K, λ). (16)
hθ(·) is any machine learning model parameterized by learnable parameters θ, such as multilayer perceptrons (MLPs). Xfea ∈ Rn×d denotes the initial node features. The model hθ(·) will first transform the initial node features as Xin = hθ(Xfea). AMP takes hθ(Xfea) as input, and performs K steps of AMP with the hyperparameter λ. Similar to the majority of existing GNN models, the training objective is the cross-entropy classification loss on the labeled nodes, and the whole model is trained in an end-to-end way. Note that AirGNN is very efficient as explained in Remark 1, and it only requires two hyperparameters K and λ without introducing additional parameters to learn, which could reduce the risk of overfitting.

heading: 4 Experiment

text: In this section, we aim to verify the effective of the proposed adaptive message passing scheme (AMP) and the AirGNN model through the semi-supervised node classification tasks. Specifically, we try to answer the following questions: (1) How does AirGNN perform on abnormal and normal features? (Section 4.2 and 4.3) and (2) How does AirGNN work by adjusting the adaptive residual? (Section 4.4)

heading: 4.1 Experimental Settings

text:
Datasets and baselines. We conduct experiments on 8 real-world datasets including three citation graphs, i.e., Cora, Citeseer, Pubmed [21], two co-authorship graphs, i.e., Coauthor CS and Coauthor Physics [22], two co-purchase graphs, i.e., Amazon Computers and Amazon Photo [22], and one OGB dataset, i.e., ogbn-arxiv [23]. Due to the space limit, we only present the results on Cora, Citeseer, and Pubmed in this section, but defer the results on other datasets to Appendix D.1. More details about the data statistics and data splits are summarized in Appendix B. The proposed AirGNN is compared with representative GNNs, including GCN [3], GAT [6], APPNP [10] and GCNII [9]. We defer the comparison with the variants of APPNP and Robust GCN [24] to Appendix D.3 and D.4 respectively.
Parameter settings. For all baselines, we follow the best hyperparameter settings in their original papers. Additionally, we tune a best residual weight α for APPNP and GCNII in the range [0, 1]. For AirGNN, we use a two-layer MLP as the base model hθ(·), following APPNP. We fix the learning rate 0.01, dropout 0.8, and weight decay 0.0005. Moreover, we set γ = 12(1−λ) as suggested by Theorem 1. We choose K = 10 and tune λ in the range [0, 1]. Adam optimizer [25] is used in all experiments. We run all experiments by 10 times, and report the mean and variance.
Evaluation setting. We assess the performance of all models under two types of abnormal feature scenarios, including noisy features and adversarial features. The abnormal features are injected to randomly selected test nodes after model training. By default, all hyperparameters are tuned according to the performance on validation sets when the dataset is clean. If tuning the hyperparameter λ of AirGNN according to the validation sets after injecting abnormal features, the performance will be even better, as discussed in Appendix D.2. The performance on clean data are showed in Appendix D.5 to demonstrate that AirGNN doesn’t need to sacrifice accuracy for better robustness against abnormal features.

heading: 4.2 Performance Comparison with Noisy Features

text:
In this subsection, we consider the abnormal features in the noisy feature scenario. Specifically, we simulate the noisy features by assigning a subset of the nodes with random features sampled from a multivariate standard Gaussian distribution. Note that the selection of noise subsets has a apparent impact on the performance since some nodes are less vulnerable to abnormal features while others are more vulnerable. To reduce such variance, we report the average performance over 10 times of random selection of the noise node sets, similar to the settings in the preliminary study in Section 2. We report the node classification test accuracy on abnormal (noisy) features and normal features in Figure 5 and Figure 6, separately, under varying noisy ratio. From these figures, we can observe:
• Figure 5 shows that AirGNN significantly outperforms all baselines on all datasets in terms of the performance on noisy nodes. This verifies that AMP is able to improve the resilience to noisy features, aligning well with the design motivation.
• Figure 6 shows that AirGNN promotes the performance on normal nodes when abnormal nodes exist. This is because AMP can remove some abnormal features which are detrimental to normal nodes.

heading: 4.3 Performance Comparison with Adversarial Features

text:
In this subsection, we consider the abnormal feature scenario when the node features are maliciously attacked by the attacker to manipulate the prediction of GNNs. We use the Nettack [26] implemented in DeepRobust3 [27], a PyTorch library for adversarial attacks and defenses, to generate the adversarial features. We randomly choose 40 test nodes as the targeted nodes, and assess the performance under increasing perturbation budgets {0, 5, 10, 20, 50, 80}, where the perturbation numbers denote the number of feature dimensions that can be manipulated. The node classification accuracy on these attacked nodes are showed in Figure 7. From these figures, we can make the following observations:
3https://github.com/DSE-MSU/DeepRobust
• AirGNN is significantly more robust against adversarially attacked features than all baselines. MLP is the most vulnerable model, which demonstrates the usefulness of graph structure information in combating against abnormal node features.
• The advantages of AirGNN over the baselines become much stronger with larger perturbation budgets. This suggests that AMP can significantly improve the resilience to abnormal features.

heading: 4.4 Adaptive Residual for Abnormal & Normal Nodes

text:
To further understand and verify how AMP and AirGNN work, we investigate the adaptive score βi for each node vi. Specifically, the average adaptive scores for abnormal nodes and normal nodes in the last layer of AMP are computed separately. In the noisy feature scenario, we fix ratio of noisy nodes as 10%. In the adversarial feature scenario, we choose 40 target nodes and fix the perturbation number as 80. The results in noisy and adversarial feature scenarios are showed in Table 1 and Table 2, respectively. From these tables, we can observe:
• On the one hand, it can be clearly observed that in both scenarios, the average adaptive scores for abnormal nodes are significantly higher than those for normal nodes. Therefore, it verifies our intuition that large adaptive scores are strongly related to abnormal features.
• On the other hand, it also implies that the residual weights (i.e., 1− βi) for abnormal nodes are much lower than those of normal nodes. This perfectly aligns with our motivation to remove abnormal features by reducing their residual connections.
The study on adaptive scores verifies how the adaptive residuals in AMP and AirGNN work as designed. It corroborates that AirGNN not only tremendously boosts the resilience to abnormal features but also provides interpretable information for anomaly detection that will be useful in many security-critical scenarios since the adaptive score serves as a good indicator of abnormal nodes. Morever, it is expected that APPNP without residual will perform well on abnormal nodes but it will sacrifice the performance on normal nodes. We provide detailed comparison with APPNP w/Res and APPNP wo/Res in Appendix D.3 to show the advantages of adaptive residual of AirGNN.

heading: 5 Related Work

text:
GNNs generalize convolutional neural networks (CNN) to graph structure data through the message passing framework [1, 2, 7]. The design of message passing and GNN architectures are majorly motivated in spectral domain [3, 4] and spatial domain [5, 6, 7, 2]. Recent works have shown that the message passing in GNNs can be regarded as low-pass graph filters [12, 13]. More generally, it has been proven that message passing in many GNNs can be uniformly derived from graph signal denoising [11, 28, 29, 30]. Classic GNNs such as GCN [3] and GAT [6] achieve their best performance with shallow models, but their performance degrades when stacking more layers, which can be partially explained through oversmoothing analyses [14, 15]. Recent works propose to use residual connections or skip connections to mitigate the oversmoothing issues, and they demonstrate the potential benefits from more feature aggregations. Examples include but not limited to DeepGCNs [31], JKNet [32], GCNII [9], APPNP [10] and DeeperGNN [20]. These models use global residual connection that can not be adaptive for each node, which significantly differ from the proposed AirGNN.
Recently, there are growing interests in reducing GNNs’ vulnerability to the graph structure noise, such as Robust GCN [24], GCN-SVD [33], Pro-GNN [34], IDGL [35], ElasticGNN [36], etc. Please refer to the comprehensive surveys [37, 38] for more details. However, how to design GNNs with strong resilience to abnormal node features remains to be developed. To the best of our knowledge, AirGNN is the first GNN model that is intrinsically robust to many types of abnormal node features by design. It improves the performance in various kinds of abnormal scenarios without needing to sacrifice clean accuracy in normal settings.

heading: 6 Conclusion

text:
In this work, we discover an intrinsic tension between feature aggregation and residual connection in the message passing scheme of GNNs, as well as the corresponding performance tradeoff between nodes with abnormal and normal features. We analyze possible reasons to explain these findings from the perspective of graph Laplacian smoothing. Our understandings further motivate us to propose a simple, efficient, interpretable and adaptive message passing scheme as well as a new GNN model with adaptive residual, named AirGNN. AirGNN provides a node-wise adaptive transition between feature aggregation and residual connection, and the significant advantages of AirGNN are demonstrated through extensive experiments.
Acknowledgments and Disclosure of Funding
This research is supported by the National Science Foundation (NSF) under grant numbers IIS1714741, CNS1815636, IIS1845081, IIS1907704, DRL2025244, IIS1928278, IIS1955285,
IOS2107215, IOS2035472 and Army Research Office (ARO) under grant number W911NF-211-0198.
Societal Impact and Limitations
The methodology proposed in this paper might have significant positive societal impact since it reduces machine learning models’ vulnerability to abnormal datasets that generally exist in real-world applications, especially in many security-critical scenarios. In practice, for a given graph, we do not have the prior knowledge about if the graph is clean, has noisy features or adversarial features. Therefore, algorithms like the proposed AirGNN that can work under both the clean and various abnormal feature settings are appealing. While we are unaware of any potential negative society impact, we point out two limitations of this work: (1) this paper focuses on abnormal node features and has not evaluated the performance of the proposed method when the dataset contains both abnormal node features and edges; (2) it is unclear how it performs on heterophilic graphs. It will be interesting to investigate these problems in future works.

title: Graph Neural Networks with Adaptive Residual