Copyright for this paper by its authors. Use permitted under Creative Commons
License Attribution 4.0 International (CC BY 4.0).

# Global Clipper: Enhancing Safety and Reliability of Transformer-based Object
Detection Models

Qutub Syed  Michael Paulitsch  Karthik Pattabiraman  Korbinian Hagn  Fabian
Oboril  Cornelius Buerkle  Kay-Ulrich Scholl  Gereon Hinz  Alois Knoll  Intel
Labs, Munich, Germany  University of British Columbia, Vancouver, Canada
Technical University of Munich, Munich, Germany

(2022)

###### Abstract

As transformer-based object detection models progress, their impact in
critical sectors like autonomous vehicles and aviation is expected to grow.
Soft errors causing bit flips during inference have significantly impacted DNN
performance, altering predictions. Traditional range restriction solutions for
CNNs fall short for transformers. This study introduces the Global Clipper and
Global Hybrid Clipper, effective mitigation strategies specifically designed
for transformer-based models. It significantly enhances their resilience to
soft errors and reduces faulty inferences to 0%. We also detail extensive
testing across over 64 scenarios involving two transformer models (DINO-DETR
and Lite-DETR) and two CNN models (YOLOv3 and SSD) using three datasets,
totalling approximately 3.3 million inferences, to assess model robustness
comprehensively. Moreover, the paper explores unique aspects of attention
blocks in transformers and their operational differences from CNNs.

##  1 Motivation

The adoption of Deep Neural Networks (DNNs) has significantly impacted various
sectors, including autonomous vehicles [1], aviation, healthcare [2], and
space exploration [3], where high safety and reliability are crucial. This has
spurred the growth of computer vision research communities focused on safe AI,
tackling areas such as out-of-distribution detection [4], adversarial
robustness and model interoperability [5]. A DNN-based computer vision model
processes images to classify objects and predict their bounding boxes.

Errors during inference can lead to faulty bounding boxes, significantly
altering system behaviour and underscoring the critical need for safer
hardware for model execution. DNN accelerators execute models at a high level
by constructing a computational graph that uses General matrix-to-matrix
multiplication (GEMM) [6] for extensive layer input and weight
multiplications. Key components in this process are the Multiply-accumulate
(MAC) units within the lower accelerator levels shown in fig. 1 [7]. MAC units
in DNN accelerators lack ECC protection, making them particularly vulnerable
to soft errors—a major reliability concern. Such errors, often caused by
radiation, chip ageing, manufacturing variations, or thermal issues [8, 9,
10], can alter intermediary computational values, leading to incorrect
inferences. Research shows that the soft error rate will increase with higher
transistor density, reduced feature sizes, and more cores [11, 12, 13]. For
example, a 100-core system of 16nm node may fail every 1.5 hours due to soft
errors [11], significantly affecting predictions as these propagate through
layers, as shown in fig. 1. Although soft errors do not cause permanent
damage, they can result in substantial reliability degradation.

This paper proposes a technique for mitigating soft errors in object detection
models at the application level. We simulated soft errors as bit flips using
PytorchALFI [14], an open-source tool that integrates large-scale fault
injection capabilities with PyTorch.

Figure 1: Abstract architecture of a DNN accelerator. The upper figure
illustrates potential soft errors resulting in bit flips within neurons or
weights at specific layers of the DNN model. The lower figure displays the
mean values of layers in non-faulty inferences compared to faulty inference
values when a bit-flip error is injected at the 50th layer of a transformer
model DINO-DETR.

Range restriction solutions effectively mitigate soft errors in CNN-based DNN
models by applying pre-calculated bounds at every activation layer, computed
using 20% of validation images to determine the minimum and maximum
restrictions [15]. However, current protective measures fall short against
soft errors in transformer-based vision models due to the complexity of their
architectures. Our analysis shows that existing solutions are inadequate,
necessitating significant enhancements in error mitigation strategies for
these advanced systems. Without such improvements, the robustness of
transformer-based models is compromised, highlighting the urgent need for more
sophisticated and tailored protection mechanisms. Transformer models [16, 17],
characterized by their self-attention and large linear layers, are
particularly vulnerable, as bit flip errors can cascade and significantly
alter predictions.

For example, injecting a single-bit flip error into the 50th self-attention
layer of the CoCo-trained DINO-DETR model, with 48M parameters [16], results
in faulty inference such as ghost objects as shown in fig. 2b. These errors,
which either create random high-confidence detections or erase them, can
disrupt systems dependent on these models for tasks like object tracking, as
shown in [18]. This underscores the significant impact of minor errors in
complex networks. However, applying existing range restrictions cannot
mitigate all ghost objects, as illustrated in fig. 2c.

We propose the Global Clipper and Global Hybrid Clipper range restriction
layers as a straightforward yet vital enhancement to mitigate the impacts of
soft errors in complex transformer-based models. These layers are implemented
within the activation and linear layers of self-attention blocks, crucial
points vulnerable to errors that can drastically affect network performance.
This strategy involves a nuanced balance: preserving the network’s ability to
process diverse data inputs while ensuring robustness against errors that
could lead to significant inaccuracies in outputs. By adding these range
restriction layers, Global Clipper effectively safeguards the network’s
functionality without compromising its learning capabilities, ensuring
sustained high performance even under challenging conditions. In the example,
when Global Clipper is added, all the false ghost objects created by fault
injection are removed as shown in fig. 2d.

In summary, the main contributions of this paper are as follows:

  1. 1.

We introduce the Global Clipper and Global Hybrid Clipper fault mitigation
techniques for transformer-based object detection models (section 4).

  2. 2.

We present a comprehensive study with fault injection experiments across CNN
and transformer models using three datasets, totalling 3.3 million inferences,
to analyze vulnerabilities in these vision systems. We show that the proposed
techniques are effective in reducing error rates from 6% to nearly 0% (section
5.3).

  3. 3.

We explore the unique characteristics of attention blocks in transformers
versus CNNs, providing insights into model vulnerabilities essential for
enhancing safety throughout the life-cycle of deployed transformer models
(section 5.4).

Figure 2: Visual example of faulty inferences on CoCo trained DINO-DETR model
due to bit-flips caused by the soft errors.

The following sections of the paper will first explore established methods for
addressing soft errors in transformer-based DNN models. Then, the paper will
detail the fault injection models under consideration ( as discussed in above
contributions - item 2). Next, our proposed solution, Global Clipper (item 1),
will be introduced, emphasizing its effectiveness in mitigating soft errors.
Subsequently, a thorough ablation study on implementing Global Clipper will be
conducted(item 2). Finally, an examination and comparison of the vulnerability
characteristics of Transformers and CNN across diverse datasets will be
provided (item 3).

##  2 Related Work

The reliability of safety-critical DNN models is assessed through various
metrics at the application and hardware levels, enhancing their safety and
reliability [19, 20, 15, 18]. Research has shown that DNNs are prone to soft
errors, with instances of single-bit flips leading to faulty inferences [21,
22, 23, 18]. Traditionally, protection from soft errors in hardware has
primarily involved error detection or correction codes (EDC or ECC) [24] for
memory and using residuals for computing. These mechanisms are commonly
implemented in high-end server-grade CPUs but less so in GPUs due to cost
considerations or typical relaxed application requirements. Other techniques
based on redundancy, like DMR (dual modular redundancy) and TMR (triple
modular redundancy), are also used. Despite these measures, accelerators
running DNNs may lack inherent protection. Furthermore, modular redundancy
techniques can be implemented through ensembles, as discussed by [25], which
may be used to detect and mitigate faults. However, these methods come with
substantial computational overhead. To address this issue, budding ensemble
solutions [26, 27] could be explored to reduce the computing overhead. Despite
their potential, these solutions have not yet been demonstrated for mitigating
or detecting soft errors.

Typical solutions for matrix operations in software, like algorithm-based
fault tolerance (ABFT), have been adapted for DNNs but are limited by the
overhead of checking large matrix multiplications typical in DNN applications
[28, 29]. Researchers have also developed DNN-specific solutions at the
application level based on range restriction solutions at the software level,
particularly for CNN models, and explored using activation patterns to detect
soft errors [30, 15, 29]. A small machine-learning module, reduced in
dimension, analyses these patterns to identify and reject erroneous inferences
[31, 32, 33, 34]. However, these methods face challenges related to
scalability and complexity.

The vulnerability of DNNs to soft errors, including a significant number of
CNN and few transformer-based models, is well-documented [10, 20, 18, 35, 36,
22]. However, previous studies have not extensively explored transformer
models in object detection or conducted detailed, large-scale fault injection
studies [37, 36]. This study aims to fill this gap by examining the resilience
of transformer architectures and exploring effective mitigation strategies
against soft errors in these architectures.

##  3 Fault Models

Figure 3: Integrating Global Clipper layers into transformer-based object
detection models’ self-attention blocks. *Ranger layers are recommended to be
added to activation functions, usually at ReLU layers, not SoftMax.

We consider soft errors in AI hardware accelerators, focusing on their impact
on system reliability. These errors, typically manifesting as single or
multiple-bit flips, can compromise data integrity by altering the model’s
weights and neurons, potentially skewing computations and decisions. Such
disruptions in deep neural network operations are illustrated in fig. 1,
highlighting the need for robustness strategies.

Parity or Error-Correcting Code (ECC) protects memory against soft errors,
particularly crucial for essential memory blocks due to the significant
overhead [38]. While ECC, especially SECDED (Single Error Correct, Double
Error Detect) code, can detect and correct single-bit errors, it is limited to
detecting two-bit errors without correction. This underscores the necessity
for other techniques beyond ECC, where minimizing multi-bit errors is
essential.

Our experimental setup injects faults as single or 10-bit flip errors during
inference, with isolated injections in either the neurons or weights of the
model, but not both simultaneously. This method ensures targeted and
straightforward fault analysis, with each inference undergoing a single fault
alteration. All models in our study employ 32-bit data types.

##  4 Global Clipper

Range restriction solutions [15, 30] effectively address bit flips caused by
soft errors in CNN-based models by focusing on activation layers where
convolutional layers attend to local image areas. This containment of
deviations within localized feature map areas helps prevent extensive errors.
However, these methods are less effective for transformer models, which employ
global attention mechanisms across extensive linear layers [39]. In
transformers, a bit flip can propagate errors throughout the multi-head
attention layers, significantly altering vector representations and impacting
predictions. This necessitates different mitigation strategies tailored to the
global processing nature of transformers.

We introduce a crucial enhancement to existing range restriction layers, as
illustrated in fig. 3, extending value monitoring and truncation from
activation layers to linear layers within self-attention blocks. This strategy
bolsters transformer architectures against soft-error-induced bit flips. The
Global Clipper truncates out-of-range values to a predefined interval before
deployment, operating at any activation or linear layer $\mathcal{l}$ with
bounds $B_{\text{lower}},B_{\text{upper}}$, as detailed in eq. 1. Similarly,
the Global Ranger restricts values within specified bounds, ensuring all layer
outputs adhere to expected ranges.

Figure 4: Lower and upper bounds for range restrictions, encompassing
activation layers and linear layers within the self-attention blocks of the
DINO-DETR model, are defined by the Global Clipper technique.

These layers can be seamlessly fused at the application level, ensuring
minimal overhead. Determining upper and lower bounds follows the approach
outlined in previous range restriction solutions like Ranger [15].
Specifically, these bounds are computed using 20% of the training dataset,
encompassing all activation and linear layers within transformer-based models
as shown in the fig. 4.

$L_{global\\_clipper}(x)=\begin{cases}0&\text{if}~{}x<B_{lower}~{}or~{}x>B_{upper}\\\
0&\text{if}~{}\text{Inf}\ \lor\ \text{NaN}\\\ x&otherwise\end{cases}$ (1)

We demonstrate the effectiveness of the Global Clipper solution using two
experiments that involve injecting faults into the convolution and linear
layers of existing models alongside our proposed method. Using the sample
images from validation datasets, we introduce bit flip errors at a random MSB
position in the second convolution layer of the ResNet50 backbone within the
DINO-DETR model. We then monitor the mean and variance at the 20th ReLU
activation layer and the 51st linear layer in the self-attention block
immediately following the ResNet50 encoder. This process, illustrated in
LABEL:fig:track_mean_var, allows us to evaluate the performance of the Global
Clipper and compare it with other mitigation strategies.

In the first experiment, we inject faults into the ReLU activation layer and
apply various mitigation solutions, as shown in LABEL:fig:track_mean_var_relu.
The Ranger method [15] confines data points within the extrapolated population
cluster space. In contrast, the Clipper [30], along with the proposed Global
Clipper, more tightly constrains the data to the region of non-faulty
inferences. Although all techniques maintain acceptable tolerances, Ranger and
Clipper tend to shift data points further from the original non-faulty
positions at linear layers.

In our second experiment, we inject faults into the linear layer and activate
various mitigation solutions. However, neither Ranger nor Clipper can confine
the values to their original positions due to the sensitivity of self-
attention layers to faults. Both Ranger and Clipper typically apply
restrictions only at the activation layer; however, self-attention includes a
SoftMax layer that cannot be similarly restricted without impairing the
functionality and accuracy of the block. Hence, introducing the Global
Clipper, shown in fig. 3, is essential for protecting attention blocks from
faults. With faults introduced into the 51st linear layer, the Global Clipper
successfully maintains data points closer to their original positions, as
demonstrated in LABEL:fig:track_mean_var_linear. While the definition of
Global Clipper resembles Ranger’s, refining the bounds is crucial as it
ensures that feature values remain within the sampled space in specific
layers. In both experiments mentioned above, fault injection demonstrated
using sample data doesn’t necessarily result in faulty predictions. However,
fault locations are sampled to visualize all data points within the area,
providing a clear explanation of Global Clipper. Additionally, more faults
affect the layer values, causing them to shift into log space compared to the
fault-free visualizations. There are a few cases in which Global Clipper may
not function out of the box on certain models, requiring slight modifications.
These exceptions and adaptations are further elaborated upon in section 5.3.

##  5 Experiments

###  5.1 Experimental Setup

As introduced in section 1, soft errors, characterized by transient bit flip
errors at the application level, impact individual inferences and last only
until the next data fetch from memory. Bit flips at the sign and most
significant bits of the mantissa minimally affect value ranges; however, flips
at the sign and exponent bits of IEEE 754 floating-point arithmetic can alter
predictions, which does not significantly change when testing other formats
like BFloat16 as seen in [21, 18].

Our study assesses vulnerabilities in two transformer-based models, DINO-DETR
[16] and Lite-DETR [17], and two CNN models, YOLOv3 [40] and SSD [41], across
the CoCo [42], KITTI [43], and BDD100K datasets [44]. We conducted over 64
experiments, totalling about 3.3 million inferences. Each model undergoes
experiments with random and targeted fault injections across all linear layers
of self-attention blocks for transformers and convolution layers for CNN
models.

Each data point extracted from these experiments includes 50,000 inferences
from 1,000 image samples with random faults and 10,000 inferences per targeted
fault experiment at each layer. Each experiment set repeats with baseline and
proposed mitigation techniques, like Clipper, Global Clipper or Global Hybrid
Clipper, allowing a thorough analysis of model vulnerabilities and the
effectiveness of mitigation strategies.

###  5.2 Evaluation Metrics

Accuracy metrics like AP50 or mAP [42] and their variants [45] are standard
for evaluating fault injections in object detection models. Occasionally,
these faults create ghost objects with lower confidence scores, not affecting
the overall AP50 due to their exclusion in the PR curve area under curve (AUC)
calculations [18]. To address this issue, the $IVMOD$ metric [18], which is
insensitive to PR curve averaging, is employed. Faults that do not alter the
model’s outcome are considered benign. In contrast, significant faults are
categorized into SDC (silent data corruption) and DUE (detectable and
unrecoverable error) as recognized by the safety and reliability community
[8].

By defining SDC and DUE [8] as critical faults, we enhance our ability to
assess vulnerabilities effectively. This study introduces these conditions
into the $IVMOD_{fd}$ metric for faulty detections (see eq. 2, eq. 3, and eq.
4). Unlike the Global Clipper, the Global Ranger restricts values without
truncating them. Vulnerability is evaluated by monitoring AP50 accuracy and
$IVMOD_{fd}$. For example, if 30 out of 100 sampled images show detection
discrepancies or encounter $NaN$ or $inf$ errors due to bit flips, the
$IVMOD_{fd}$ would be 30%. Additionally, $IVMOD_{fd}$ can be used to estimate
DNN accelerator vulnerability in terms of FIT rates [8] and other risk factors
throughout the hardware’s lifecycle [19, 20]. However, these aspects are
beyond this paper’s scope. Hereafter, $IVMOD_{fd}$ will be interchangeably
called faulty detections. Moreover, in this study, $IVMOD_{fd}$ solely
considers the 9 higher-order bits, including the sign and exponent bits, as
described in section 3.

$\begin{split}IVMOD_{SDC}=\frac{1}{N}\sum^{N}_{i=1}[(FP_{\text{orig}})_{i}\neq(FP_{\text{corr}})_{i}\lor\
\\\ (FN_{\text{orig}})_{i}\neq(FN_{\text{corr}})_{i}\big{]}\\\ \end{split}$
(2)

$\begin{split}IVMOD_{DUE}=\frac{1}{N}\sum^{N}_{i=1}[\text{Inf}\ \lor\
\text{NaN}]\\\ \end{split}$ (3)

$\begin{split}IVMOD_{fd}=IVMOD_{SDC}\lor\ IVMOD_{DUE},\\\ \end{split}$ (4)

###  5.3 Results: Global Clipper on Transformer models

This section showcases the superior mitigation capability of Global Clipper
and Global Hybrid Clipper over existing solutions like Ranger and Clipper.
This is demonstrated on DINO-DETR and Lite-DETR models with various datasets
injected with single bit-flip errors (see LABEL:fig:dino_detr_sdc_ap50 and
LABEL:fig:lite_detr_sdc_ap50), we also investigate vulnerability across
datasets, including layer-wise fault injection experiments.

The LABEL:fig:dino_detr_sdc_ap50 illustrates the model’s vulnerability and
evaluates the effectiveness of the various range restriction techniques in
mitigating faults using $IVMOD_{fd}$ and AP50 metrics. This includes results
from experiments involving fault injections into weights and neurons. The
Global Clipper significantly outperforms the existing solutions like Ranger
and Clipper in mitigating the impact of the faults. The Global Clipper
performs better in mitigating faults occurring in neurons, reducing faulty
detections to nearly 0%, and for weight faults, the vulnerabilities are
reduced to less than $\sim$ 3% by outperforming other state-of-the-art
algorithms like Ranger and Clipper. Additionally, as shown in
LABEL:fig:dino_detr_ap50, a single bit flip in the inferences notably impacts
the AP50 metric. For instance, the AP50 of DINO-DETR trained on KITTI
decreases from 86.9 to 81.6 when injected with weight faults, and the Global
Clipper effectively restores the AP50 to its original accuracy. In the context
of DINO-DETR, Global Ranger and Global Clipper demonstrate similar mitigation
performance. However, their performance is not as expected when applied to
Lite-DETR due to the unique transformer architecture based on DINO-DETR. This
presents an interesting scenario, leading to the introduction of the Global
Hybrid Clipper. The Global Hybrid Clipper merges Global Clipper and Ranger
layers. In this setup, Global Clipper is applied to Activation layers, while
Global Ranger is used for the linear layers within self-attention blocks. In
this scenario, the hybrid version of Global Clipper restores performance
accuracy to baseline and reduces faulty detections significantly, as shown in
LABEL:fig:lite_detr_sdc_ap50. For instance, it decreases from 4.5% to 0.5% in
the case of neuron faults.

Figure 6: Individual components of Encoders and Decoders (Attention Blocks)

###  5.4 Ablation Study

####  5.4.1 Vulnerability and Resiliency Analysis of Attention and Convolution
Layers

This section outlines the ablation study, including single-bit and multiple-
bit (10-bit) flip experiments on CNN and transformer-based object detection
models. Results are shown in LABEL:fig:single_bit_flip_sdc and
LABEL:fig:multi_bit_flip_sdc. Global Clipper and Global Hybrid Clipper were
employed for transformers, while Clipper was used for CNNs due to its superior
performance over Ranger (see LABEL:fig:CNN_1FI and LABEL:fig:CNN_10FI).

Each plot represents 10,000 inferences, highlighting Global Clipper’s superior
handling of bit-flip errors. Comparing CNN and self-attention layers reveals
that transformers generally exhibit greater fault-injection resilience than
CNNs. For instance, DINO-DETR and Lite-DETR consist of six encoders and
decoders, each with four linear layers (as depicted in fig. 6), resulting in
48 layers per model. These transformer layers exhibit distinct characteristics
in response to fault injections compared to CNNs. These observations are
consistent across both single-bit and multi-bit flip experiments
(LABEL:fig:single_bit_flip_sdc and LABEL:fig:multi_bit_flip_sdc). Depending on
the transformer model variant, encoders and decoders are interchangeable with
self-attention blocks.

CNNs show no discernible pattern across layers, suggesting uniform
susceptibility to generating faulty detections under bit-flip errors. In
contrast, transformer models display greater inherent resilience. With
suitable mitigation techniques, transformers offer enhanced safety against
soft errors compared to CNN-based object detection models, making them
preferable for applications demanding robustness against bit-flip errors.

Figure 8: Integrating Global Clipper in Attention Block with minimal overhead.

To better understand the transformer’s vulnerability, we analyzed single-bit
flip errors in attention block linear layers A, B, C, and D (see fig. 6). Data
from LABEL:fig:transform_1FI and LABEL:fig:transform_10FI was segmented into
four plots in LABEL:fig:attention_block_vulnerability, illustrating each
layer’s vulnerability across encoders and decoders in 12 Attention Blocks.
Layers such as Sampling Offset, Attention Weights, Value Projection, and
Output Projection showed consistent vulnerability across CoCo, KITTI, BDD100K
datasets, and DINO-DETR and Lite-DETR models for neuron faults. Variations in
vulnerability due to weight faults depend on the attention block type, as
shown in LABEL:fig:attention_block_vulnerability_DF_KDA. DINO-DETR uses
Deformable Attention (DF) Layers, while Lite-DETR employs Key-Aware Deformable
Attention (KDA) Layers, enhancing efficiency and attention mapping. The
effects of weight faults on DF and KDA blocks vary, influencing encoders and
decoders differently (LABEL:fig:attention_block_vulnerability_DF_KDA_weights),
but neuron fault vulnerability remains stable across datasets and models
(LABEL:fig:attention_block_vulnerability_DF_KDA_neurons). These observations
underscore two key transformer traits in vulnerability, distinct from CNNs:

  * •

The vulnerability estimation of a transformer model’s layers demonstrates
consistent characteristics across different datasets during inference, unlike
CNNs.

  * •

The model’s vulnerability estimation is influenced by the Self-Attention Block
variant used in the model and remains consistent across different
architectures and datasets during inference.

Thorough vulnerability analysis can greatly improve online safety and risk
management, assisting in dynamic risk assessment and model monitoring across
the lifecycle of deployed transformer-based models. Our findings demonstrate
that despite continual learning and weight adjustments, these model’s
vulnerabilities remain stable, ensuring robust and reliable AI systems are
maintained.

####  5.4.2 Integrating Global Clipper with minimal additional overhead

As depicted in fig. 3, the Global Clipper layers are integrated into four
linear layers within the Attention Block (see fig. 6). The efficiency of
adding these layers can be further evaluated by conducting a simple experiment
that examines the impact of each Global Clipper Layer. This experiment entails
the injection of a single-bit flip error at the initial stage of the Attention
Block, particularly at the Sampling Offsets layer (stage A). Following this,
the Global Clipper layer is selectively activated at different combinations of
stages A, B, C, and D within the four linear layers of the Attention Block, as
detailed in fig. 6. Our observations indicate that integrating Global Clipper
solely into the Output Projection and Value Projection layers (D and C in
Figure 6), in addition to other activation layers, yields mitigation
performance comparable to that achieved by incorporating it into all linear
layers, as illustrated in Figure 8

##  6 Conclusion

This study introduces Global Clipper and Global Hybrid Clipper to enhance the
safety of transformer-based object detection models in critical settings,
effectively minimizing faulty inferences to nearly 0%. We evaluated these
solutions by conducting fault injection campaigns with transformer and CNN
models across three datasets, totalling approximately 3.3 million inferences.
Our extensive experiments and findings indicate that transformer models
exhibit better inherent resilience to soft errors than CNN models. Evaluating
these solutions provides insights into their effectiveness in real-world
applications, contributing significantly to model robustness and computer
vision safety. Future research should explore these solutions in transformer-
based semantic segmentation and video tracking models to further enhance
safety.

## References

  * Wang et al. [2020] J. Wang, L. Zhang, Y. Huang, J. Zhao, F. Bella,  Safety of autonomous vehicles,  Journal of advanced transportation 2020 (2020) 1–13.
  * Habli et al. [2020] I. Habli, T. Lawton, Z. Porter,  Artificial intelligence in health care: accountability and safety,  Bulletin of the World Health Organization 98 (2020) 251.
  * Oche et al. [2021] P. A. Oche, G. A. Ewa, N. Ibekwe,  Applications and challenges of artificial intelligence in space missions,  IEEE Access (2021).
  * Gawlikowski et al. [2023] J. Gawlikowski, C. R. N. Tassi, M. Ali, J. Lee, M. Humt, J. Feng, A. Kruspe, R. Triebel, P. Jung, R. Roscher, et al.,  A survey of uncertainty in deep neural networks,  Artificial Intelligence Review 56 (2023) 1513–1589.
  * Xu et al. [2019] F. Xu, H. Uszkoreit, Y. Du, W. Fan, D. Zhao, J. Zhu,  Explainable ai: A brief survey on history, research areas, approaches and challenges,  in: Natural Language Processing and Chinese Computing: 8th CCF International Conference, NLPCC 2019, Dunhuang, China, October 9–14, 2019, Proceedings, Part II 8, Springer, 2019, pp. 563–574.
  * Gao et al. [2023] J. Gao, W. Ji, F. Chang, S. Han, B. Wei, Z. Liu, Y. Wang,  A systematic survey of general sparse matrix-matrix multiplication,  ACM Computing Surveys 55 (2023) 1–36.
  * Chen et al. [2020] Y. Chen, Y. Xie, L. Song, F. Chen, T. Tang,  A survey of accelerator architectures for deep neural networks,  Engineering 6 (2020) 264–274.
  * Mukherjee et al. [2005] S. S. Mukherjee, J. Emer, S. K. Reinhardt,  The soft error problem: An architectural perspective,  in: 11th International Symposium on High-Performance Computer Architecture, IEEE, 2005, pp. 243–247.
  * Haque and Pande [2010] I. S. Haque, V. S. Pande,  Hard data on soft errors: A large-scale assessment of real-world error rates in gpgpu,  in: 2010 10th IEEE/ACM International Conference on Cluster, Cloud and Grid Computing, 2010, pp. 691–696. doi:10.1109/CCGRID.2010.84.
  * Ibrahim et al. [2020] Y. Ibrahim, H. Wang, J. Liu, J. Wei, L. Chen, P. Rech, K. Adam, G. Guo,  Soft errors in dnn accelerators: A comprehensive review,  Microelectronics Reliability 115 (2020) 113969.
  * Upasani [2016] G. R. Upasani,  Soft error mitigation techniques for future chip multiprocessors (2016).
  * Borkar [2005] S. Borkar,  Designing reliable systems from unreliable components: the challenges of transistor variability and degradation,  Ieee Micro 25 (2005) 10–16.
  * Shivakumar et al. [2002] P. Shivakumar, M. Kistler, S. W. Keckler, D. Burger, L. Alvisi,  Modeling the effect of technology trends on the soft error rate of combinational logic,  in: Proceedings International Conference on Dependable Systems and Networks, IEEE, 2002, pp. 389–398.
  * Gräfe et al. [2023] R. Gräfe, Q. S. Sha, F. Geissler, M. Paulitsch,  Large-scale application of fault injection into pytorch models-an extension to pytorchfi for validation efficiency,  in: 2023 53rd Annual IEEE/IFIP International Conference on Dependable Systems and Networks-Supplemental Volume (DSN-S), IEEE, 2023, pp. 56–62.
  * Chen et al. [2021] Z. Chen, G. Li, K. Pattabiraman,  A low-cost fault corrector for deep neural networks through range restriction,  in: 2021 51st Annual IEEE/IFIP International Conference on Dependable Systems and Networks (DSN), IEEE, 2021, pp. 1–13.
  * Zhang et al. [2022] H. Zhang, F. Li, S. Liu, L. Zhang, H. Su, J. Zhu, L. M. Ni, H.-Y. Shum,  Dino: Detr with improved denoising anchor boxes for end-to-end object detection,  arXiv preprint arXiv:2203.03605 (2022).
  * Li et al. [2023] F. Li, A. Zeng, S. Liu, H. Zhang, H. Li, L. Zhang, L. M. Ni,  Lite detr: An interleaved multi-scale encoder for efficient detr,  in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2023, pp. 18558–18567.
  * Qutub et al. [2022] S. Qutub, F. Geissler, Y. Peng, R. Gräfe, M. Paulitsch, G. Hinz, A. Knoll,  Hardware faults that matter: Understanding and estimating the safety impact of hardware faults on object detection dnns,  in: International Conference on Computer Safety, Reliability, and Security, Springer, 2022, pp. 298–318.
  * Neale and Sachdev [2016] A. Neale, M. Sachdev,  Neutron radiation induced soft error rates for an adjacent-ecc protected sram in 28 nm cmos,  IEEE Transactions on Nuclear Science 63 (2016) 1912–1917.
  * Li et al. [2017] G. Li, S. K. S. Hari, M. Sullivan, T. Tsai, K. Pattabiraman, J. Emer, S. W. Keckler,  Understanding error propagation in deep learning neural network (dnn) accelerators and applications,  in: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2017, pp. 1–12.
  * Geissler et al. [2021] F. Geissler, S. Qutub, S. Roychowdhury, A. Asgari, Y. Peng, A. Dhamasia, R. Graefe, K. Pattabiraman, M. Paulitsch,  Towards a safety case for hardware fault tolerance in convolutional neural networks using activation range supervision,  arXiv preprint arXiv:2108.07019 (2021).
  * Geissler et al. [2023] F. Geissler, S. Qutub, M. Paulitsch, K. Pattabiraman,  A low-cost strategic monitoring approach for scalable and interpretable error detection in deep neural networks,  in: International Conference on Computer Safety, Reliability, and Security, Springer, 2023, pp. 75–88.
  * Asgari Khoshouyeh et al. [2023] A. Asgari Khoshouyeh, F. Geissler, S. Qutub, M. Paulitsch, P. Nair, K. Pattabiraman,  Structural coding: A low-cost scheme to protect cnns from large-granularity memory faults,  in: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2023, pp. 1–17.
  * Hamming [1950] R. W. Hamming,  Error detecting and error correcting codes,  The Bell system technical journal 29 (1950) 147–160.
  * Shlezinger et al. [2021] N. Shlezinger, E. Farhan, H. Morgenstern, Y. C. Eldar,  Collaborative inference via ensembles on the edge,  in: ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2021, pp. 8478–8482.
  * Qutub et al. [2023] S. S. Qutub, N. K. Cihangir, R. Rosales, M. Paulitsch, K. Hagn, F. Geissler, Y. Peng, G. Hinz, A. C. Knoll,  Bea: Revisiting anchor-based object detection dnn using budding ensemble architecture.,  in: BMVC, 2023, pp. 792–797.
  * Qutub et al. [2024] S. S. Qutub, M. Paulitsch, K.-U. Scholl, N. K. Cihangir, K. Hagn, F. Oboril, G. Hinz, A. Knoll,  Situation monitor: Diversity-driven zero-shot out-of-distribution detection using budding ensemble architecture for object detection,  in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) Workshops, 2024, pp. 3502–3511.
  * Chen [2013] Z. Chen,  Online-abft: An online algorithm based fault tolerance scheme for soft error detection in iterative methods,  ACM SIGPLAN Notices 48 (2013) 167–176.
  * Zhao et al. [2020] K. Zhao, S. Di, S. Li, X. Liang, Y. Zhai, J. Chen, K. Ouyang, F. Cappello, Z. Chen,  Ft-cnn: Algorithm-based fault tolerance for convolutional neural networks,  IEEE Transactions on Parallel and Distributed Systems 32 (2020) 1677–1689.
  * Hoang et al. [2020] L.-H. Hoang, M. A. Hanif, M. Shafique,  Ft-clipact: Resilience analysis of deep neural networks and improving their fault tolerance using clipped activation,  in: 2020 Design, Automation & Test in Europe Conference & Exhibition (DATE), IEEE, 2020, pp. 1241–1246.
  * Schorn and Gauerhof [2020] C. Schorn, L. Gauerhof,  Facer: A universal framework for detecting anomalous operation of deep neural networks,  in: 2020 IEEE 23rd International Conference on Intelligent Transportation Systems (ITSC), IEEE, 2020, pp. 1–6.
  * Schorn et al. [2018] C. Schorn, A. Guntoro, G. Ascheid,  Efficient on-line error detection and mitigation for deep neural network accelerators,  in: Computer Safety, Reliability, and Security: 37th International Conference, SAFECOMP 2018, Västerås, Sweden, September 19-21, 2018, Proceedings 37, Springer, 2018, pp. 205–219.
  * Zhao et al. [2022] F. Zhao, C. Zhang, N. Dong, Z. You, Z. Wu,  A uniform framework for anomaly detection in deep neural networks,  Neural Processing Letters 54 (2022) 3467–3488.
  * Mahmoud et al. [2020] A. Mahmoud, S. K. S. Hari, C. W. Fletcher, S. V. Adve, C. Sakr, N. Shanbhag, P. Molchanov, M. B. Sullivan, T. Tsai, S. W. Keckler,  Hardnn: Feature map vulnerability evaluation in cnns,  arXiv preprint arXiv:2002.09786 (2020).
  * Xue et al. [2023] X. Xue, C. Liu, Y. Wang, B. Yang, T. Luo, L. Zhang, H. Li, X. Li,  Soft error reliability analysis of vision transformers,  IEEE Transactions on Very Large Scale Integration (VLSI) Systems (2023).
  * Roquet et al. [2024] L. Roquet, F. F. dos Santos, P. Rech, M. Traiola, O. Sentieys, A. Kritikakou,  Cross-layer reliability evaluation and efficient hardening of large vision transformers models (2024).
  * Agarwal et al. [2023] U. K. Agarwal, A. Chan, K. Pattabiraman,  Resilience assessment of large language models under transient hardware faults,  in: 2023 IEEE 34th International Symposium on Software Reliability Engineering (ISSRE), IEEE, 2023, pp. 659–670.
  * Lotfi et al. [2019] A. Lotfi, S. Hukerikar, K. Balasubramanian, P. Racunas, N. Saxena, R. Bramley, Y. Huang,  Resiliency of automotive object detection networks on GPU architectures,  Proceedings - International Test Conference 2019-Novem (2019) 1–9. doi:10.1109/ITC44170.2019.9000150.
  * Dosovitskiy et al. [2020] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, et al.,  An image is worth 16x16 words: Transformers for image recognition at scale,  arXiv preprint arXiv:2010.11929 (2020).
  * Redmon and Farhadi [2018] J. Redmon, A. Farhadi,  Yolov3: An incremental improvement,  arXiv preprint arXiv:1804.02767 (2018).
  * Liu et al. [2016] W. Liu, D. Anguelov, D. Erhan, C. Szegedy, S. Reed, C.-Y. Fu, A. C. Berg,  Ssd: Single shot multibox detector,  in: Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11–14, 2016, Proceedings, Part I 14, Springer, 2016, pp. 21–37.
  * Lin et al. [2014] T.-Y. Lin, M. Maire, S. Belongie, J. Hays, P. Perona, D. Ramanan, P. Dollár, C. L. Zitnick,  Microsoft coco: Common objects in context,  in: Computer Vision–ECCV 2014: 13th European Conference, Zurich, Switzerland, September 6-12, 2014, Proceedings, Part V 13, Springer, 2014, pp. 740–755.
  * Geiger et al. [2012] A. Geiger, P. Lenz, R. Urtasun,  Are we ready for autonomous driving? the kitti vision benchmark suite,  in: 2012 IEEE conference on computer vision and pattern recognition, IEEE, 2012, pp. 3354–3361.
  * Yu et al. [2020] F. Yu, H. Chen, X. Wang, W. Xian, Y. Chen, F. Liu, V. Madhavan, T. Darrell,  Bdd100k: A diverse driving dataset for heterogeneous multitask learning,  in: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 2020, pp. 2636–2645.
  * dos Santos et al. [2019] F. F. dos Santos, P. Navaux, L. Carro, P. Rech,  Impact of reduced precision in the reliability of deep neural networks for object detection,  in: 2019 IEEE European Test Symposium (ETS), IEEE, 2019, pp. 1–6.


# Urban hierarchy and spatial diffusion over the innovation life cycle

Eszter Bokányi1,2,∗, Martin Novák1,3, Ákos Jakobi4,5, Balázs Lengyel1,2

(1ELKH Centre for Economic and Regional Studies; Agglomeration and Social
Networks Lendület Research Group,
Budapest, H-1097, Hungary
2Corvinus University of Budapest; Laboratory for Networks, Technology and
Innovation; Budapest, H-1093, Hungary
3ELKH Wigner Research Centre for Physics; Budapest, H-1121, Hungary
4Eötvös Loránd University; Department for Regional Science; Budapest, H-1117,
Hungary
5ELKH Research Centre for Astronomy and Earth Sciences, Spatial Big Data Lab;
Budapest, H-1121, Hungary
∗Corresponding author<EMAIL_ADDRESS>
)

###### Abstract

Successful innovations achieve large geographical coverage by spreading across
settlements and distances. For decades, spatial diffusion has been argued to
take place along the urban hierarchy such that the innovation first spreads
from large to medium cities then later from medium to small cities. Yet, the
role of geographical distance, the other major factor of spatial diffusion,
was difficult to identify in hierarchical diffusion due to missing data on
spreading events. In this paper, we exploit spatial patterns of individual
invitations on a social media platform sent from registered users to new users
over the entire life cycle of the platform. This enables us to disentangle the
role of urban hierarchy and the role of distance by observing the source and
target locations of flows over an unprecedented timescale. We demonstrate that
hierarchical diffusion greatly overlaps with diffusion to close distances and
these factors co-evolve over the life cycle; thus, their joint analysis is
necessary. Then, a regression framework is applied to estimate the number of
invitations sent between pairs of towns by years in the life cycle with the
population sizes of the source and target towns, their combinations, and the
distance between them. We confirm that hierarchical diffusion prevails
initially across large towns only but emerges in the full spectrum of
settlements in the middle of the life cycle when adoption accelerates. Unlike
in previous gravity estimations, we find that after an intensifying role of
distance in the middle of the life cycle a surprisingly weak distance effect
characterizes the last years of diffusion. Our results stress the dominance of
urban hierarchy in spatial diffusion and inform future predictions of
innovation adoption at local scales.

##  1 Introduction

Fuelled by the recognition that innovation drives economic progress [1, 2],
the adoption dynamics of new technologies and products over their life cycle
have been analyzed and modeled from the middle of the $20^{\mathrm{th}}$
century [3, 4]. Griliches [5] quickly made geography an important field of
this discussion by demonstrating that adoption is faster in the proximity of
the original location. Later, Hägerstand [6] modeled spatial diffusion with a
special focus on cities; he distinguished categories of hierarchical (from a
larger to a smaller city) and neighborhood (determined by geographical
proximity) processes. This classic work was applied in a wide range of
contexts including the diffusion of ideas [7], fashion [8], culture [9], or
e-shopping [10]. The mixture of hierarchical and neighborhood diffusion
mechanisms are thought to resemble routing behavior in social media [11] such
that spreading initially occurs over great distances between large cities but
becomes more and more local when it reaches smaller towns. However, the
separation of hierarchical versus neighborhood diffusion is still problematic,
because large cities are typically surrounded by smaller towns [12], which
demands further analysis with data on spreading events.

More recently, adoption dynamics in social networks have received growing
attention [13, 14]. In this stream of literature, diffusion is modeled as a
complex contagion process, in which the probability of individuals’ adoption
increases as network neighbors adopt [15]. This model can predict the online
spreading of behavior and the diffusion of online products at the scale of
countries well [16, 17, 18]. However, local predictions of the complex
contagion framework still suffer from systemic biases [19]. City size and
distance from the origin have remained major confounding factors in modeling
spatial adoption [20] and still limit marketing applications at local scales
[21].

In this paper, we aim to better understand the spatial diffusion of
innovations by separating the role of urban hierarchy and geographical
distance. This is done by analyzing a unique database retrieved from a social
media platform called iWiW, the most popular online social network in Hungary
in the pre-Facebook era, adopted by 30% of the total population. We can trace
the platform’s spatial diffusion through the dynamics of geolocated platform
registrations, and the data also allows us to observe each spreading event -
that are accepted invitations registered users sent to new users - over the
whole life cycle of the product from 2002 to 2012. Thus, we can reconstruct
the direction of flows with unprecedented spatial granularity, scale, and time
horizon.

We demonstrate the problems of hierarchical versus neighborhood diffusion
categorization by showing their overlaps. In the middle of the life cycle,
when adoption accelerates, the most probable spreading path across towns form
long chains signaling that hierarchical diffusion emerges across the full
spectrum of towns. In this period, large size differences between source and
target towns of spreading events coincide with large shares of invitations
sent to very small distances in agglomeration areas where invitations tend to
cascade from large cities to neighboring small towns.

To disentangle hierarchical and neighborhood diffusion, we estimate the number
of invitations sent between pairs of towns by years with zero-inflated
negative binomial regressions, in which the population sizes of the source and
the target towns, their combination, and the distance between them are the
explanatory variables. This framework allows us to separate the single effect
of city size and distance and also consider their joint effect. By considering
the direction of flows, we discover that hierarchical diffusion prevails
initially across large towns only but emerges in the full spectrum of
settlements in the middle of the life cycle. Interestingly, we discover a
characteristic size difference as a sign of hierarchical diffusion at the end
of the life cycle. Fixing the urban hierarchy variables at given levels, we
find that the role of distance is minor in the beginning, intensifies in the
middle, and shades away towards the end of the life cycle. This finding
contradicts previous gravity estimations that have suggested monotonous
intensification of distance effect.

Our results confirm that spatial diffusion of innovations occurs along the
hierarchical order of cities. We provide new insights into the dynamics of the
exact mechanism over the innovation life cycle. This new evidence is important
for local adoption predictions that aim to overcome the systemic biases of
settlement size and geographical distance in network diffusion models.

##  2 Results

###  2.1 Data

Our analysis is based on a dataset collected from iWiW, a Hungarian online
social networking site launched in 2002. The service was a remarkable
innovation in the Hungarian market of online services of its time and became
highly popular in 2005. Almost 3 million profiles were registered until 2008
(30% of total population) when iWiW was the most visited Hungarian webpage.
The life cycle of iWiW is depicted by monthly registrations in Supporting
Information 1, in which registrations are broken down to age groups as well.
The website lost its primacy from 2010, when Facebook overtook the market and
the iWiW owner stopped the service in 2014. We have access to the anonymized
version of all public profile data on iWiW, including friendship ties,
demographics, location, registration date, and date of last login. Data was
collected and delivered by the data owner in January 2013. Previous research
on iWiW has looked into the spatial pattern of the iWiW social network and its
relation with socio-economic phenomena [22, 23, 24], churn from the website
[25, 26], and diffusion [20, 27].

In this paper, we analyze invitation patterns that capture directed and
geolocated spreading events at unprecedented scale and time horizon. Until
2012, users could only register on the website if they have been invited by
another iWiW user. The individual-level data used in this paper contains the
ID of iWiW users, the ID of those users who sent the invitation, self-reported
location information of users at the settlement level, and day of
registration. Thus, we are able to construct a network of invitations from our
dataset, where nodes are the iWiW users, and directed edges between them
represent the invitor-invited relationships at the individual level. This way,
we create a network with 3,059,363 users and 2,636,779 invitations between
them. We then aggregate our individual-level network into the 2555 Hungarian
settlements given in the users town of residence field in each year of the
life cycle. In this dynamic settlement-level spatial network, we have directed
edges with an edge weight corresponding to the number of invitations going
from the source settlement to the target settlement in the given year. We omit
self-loops from the analysis, where the settlement of the source node and the
target node are identical.

###  2.2 Overlap of hierarchical and neighborhood diffusion

Figure 1: Hierarchical and neighborhood diffusion over the innovation life
cycle. (a) Most probable invitation paths given by the solution of the minimum
weight branching problem on the transformed network with $\tilde{p}_{ST}$
weights. Nodes correspond to settlements. (b) Minimum weight invitation trees
of the first row with settlements positioned on the map of Hungary. Edges are
colored according to the size of the source settlement (edges with Budapest as
a source are blue), that is also indicated by the size of the nodes. (c)
Distribution of tree edges with respect to two measures: distance between
source and target settlements (horizontal axis), and $\log_{10}$ of source and
target population size fraction (vertical axis). The vertical black line
separates tree edges with less than 20 km distance between source and target,
the horizontal black line separates tree edges that go down the settlement
hierarchy (source size is larger than the target size), and that go up the
settlement hierarchy (target size is larger than the source size). The top
left quadrant corresponds to edges that are both sent to a very close distance
and in a downwards hierarchical pattern.

To demonstrate the problem that hierarchical and neighborhood diffusion often
overlap, Figure 1 illustrates the most probable invitation tree in the spatial
diffusion measured from the settlement-level network across the iWiW life
cycle. We denote the number of invitations in a given year going from
settlement $S$ (source) to settlement $T$ (target) by $w_{ST}$. To identify
most probable diffusion paths, we normalize $w_{ST}$ by the number of total
outgoing invitations such that edges represent outgoing invitation
probabilities $p_{ST}$, where $p_{ST}=\frac{w_{ST}}{\sum_{T}w_{ST}}$. The
probability of a cross-settlement diffusion path is the product of these
$p_{ST}$ values along the path edges. We normalize weights by taking the
logarithm and multiply by -1 so that the new weights are
$\tilde{p}_{ST}=-\log\frac{w_{ST}}{\sum_{T}w_{ST}}$. This transformation makes
the maximum probability tree of the graph the one for which the sum of the
weights is minimal. Using the modified weights enables us to search for the
most probable invitation tree [28, 29] by detecting the minimum weight
branching in the invitation network [30]. We use the efficient implementation
of the Chu-Liu-Edmonds algorithm for obtaining the minimum weight spanning
tree by Gabow [31].

Figure 1a reports a characteristic change in the pattern of major flows of
diffusion across settlements. In 2003-2005, the paths in the tree are short
such that leaves can be reached through only a low number of internal
vertices. The tree even includes a node that has an especially large degree
corresponding to Budapest that acts as a major source of invitations during
this period. These observations suggest fast diffusion cascades, in which few
settlements stand out in distributing the innovation. However, the pattern
changes in 2006-2008 when paths in the diffusion tree form long chains, and
large-degree nodes disappear signaling a smoother and more equal diffusion
through the whole size-spectrum of settlements. This tree with long paths
remains in years 2009-2012 as well but smaller cascades emerge again at the
tip of the branches.

The spatial pattern of these diffusion trees illustrated in Figure 1b,
provides support to the hierarchical diffusion of Hägerstrand [6].
Hierarchical diffusion of iWiW initially happened only from the largest cities
that have mediated the innovation to all other settlements even to large
distances. Later, from the middle of the life cycle, the long diffusion paths
have emerged along the urban hierarchy such that flows of spreading went from
slightly larger to slightly smaller settlements to a smaller geographical
distance. Finally, at the end of the life cycle, both mechanisms contribute to
the diffusion.

Next, we further demonstrate the co-evolution of the hierarchical and
neighborhood diffusion over the life cycle in Figure 1c. This is done by
plotting the geographical distance between settlements connected by the
diffusion tree in Figure 1a against the population difference (measured as the
$\log_{10}$ of the fraction of the source and the target population sizes).
The first period is characterized by long-distance cascades of hierarchical
diffusion. This changes to local diffusion across similar settlements in the
second period. Finally, the emerging small cascades remain local.

Taken together, in the above descriptive analysis we find by looking at the
most probable diffusion tree across settlements that hierarchical diffusion
and neighborhood diffusion are closely related and co-evolve over the life
cycle. Thus, their joint systematic analysis is necessary. Supplementary
Information 2 further describes the dynamics of hierarchical and neighborhood
invitations.

###  2.3 Disentangling hierarchical and neighborhood diffusion

According the observations presented above, the effect of urban hierarchy in
spatial diffusion must be disentangled from other major geographical factors.
The first of such geographical factors is the effect of the urban
concentration captured by urban scaling law and the second is the effect of
distance captured by the gravity law, both of which have been previously found
to characterize spatial diffusion of iWiW [20]. The scaling effect [32]
manifests in disproportionally more invitations sent from large settlements
early on in the life cycle compared to what their population size would
linearly suggest, but this superlinear scaling disappears throughout the
majority and the laggard phases [20]. The gravity effect captures the
intensity of diffusion between settlements $S$ and $T$ as a function of
populations $P_{S}$ and $P_{T}$ and the term $d^{-\chi}$, where $d$ is the
Euclidean distance between the two settlements, and $\chi$ is the gravity
coefficient, which is found to have an increasing power on diffusion over the
life cycle [20].

We apply a regression framework to capture the role of hierarchical diffusion,
the effect that an innovation is more likely to go from larger settlements
towards smaller settlements, assuming the co-presence of urban scaling of
adoption and the gravity law in sending invitations. In every year, we observe
the number of invitations going between settlement pairs $w_{ST}$, but there
are even more settlement pairs between which there is no invitation process at
all. To account for these excess zero observations, we use a zero-inflated
negative binomial regression model [33]. In this regression technique, we
assume that the excess zero counts come from a logit model and the remaining
counts come from a negative binomial model specified in Materials and Methods.
We can predict invitation counts by including all possible settlement pairs
into the data.

We estimate the following equation for the invitation counts $w_{ST}$:

$\log w_{ST}\approx\alpha(\log P_{S})^{2}+\beta\cdot(\log P_{S}\cdot\log
P_{T})+\gamma(\log P_{T})^{2}+\delta\log P_{S}+\varepsilon\log P_{T}+\chi\log
d_{ST}+C,$ (1)

where the coefficients $\alpha$, $\beta$, $\gamma$, $\delta$, $\varepsilon$,
and $\chi$, and the constant $C$ are given by the maximum likelihood estimator
of the ZINB model in Stata [34]. To account for the gravity effect, we include
the distance term $\log d$, and the population terms corresponding to the
source and target populations $\log P_{S}$ and $\log P_{T}$. However, Figure 1
suggests that the role of hierarchy changes in the life cycle, because the
most probable tree of invitations shifts from larger towns to long hierarchy
chains, then to smaller towns dominating the invitations. Therefore, we add
second-order terms $(\log P_{S})^{2}$, $(\log P_{T})^{2}$, and $\log
P_{S}\cdot\log P_{T}$ to the estimation, which can capture and test multiple
hierarchical invitation scenarios by using terms that are able to describe a
general second-order surface in the variables $\log P_{S}$ and $\log P_{T}$.
For the full regression tables and the logit part of the model, see regression
coefficients and parameters in the tables of Supplementary Information 4.

We start with the effect of settlement size. We would like to understand how
the source and target settlement populations influence invitation counts if we
keep the distance term constant. In other words, given two settlement pairs
for which the distance is equal, how does the invitation count differ for
different $P_{S}$ and $P_{T}$ pairs. Because of the second-order terms in the
regression, the surface defined by the population terms

$\alpha(\log P_{S})^{2}+\beta\cdot(\log P_{S}\cdot\log P_{T})+\gamma(\log
P_{T})^{2}+\delta\log P_{S}+\varepsilon\log P_{T}+C$

is a general second-order surface. For us, the important range from this
surface that is defined by the regression coefficients is where both $P_{S}$
and $P_{T}$ are within realistic population sizes, that is, if they are larger
than a small village ($\approx 200$ people), but smaller than the capital city
Budapest ($\approx 2\cdot 10^{6}$ people).

Figure 2: Settlement size in spatial diffusion. The effect of source (a-c) and
target (d-f) settlement size with all other variables fixed in the three
selected years 2003, 2006, and 2011.

First, we show how $\log w_{ST}$ changes with $P_{S}$ for four different
target population sizes $P_{T}$ (Figure 2a-c) in three selected years (2003,
2006, and 2011) that are in three different phases of the life cycle. These
four different target populations roughly represent a village ($10^{3}$), a
small town ($10^{4}$), a larger town ($10^{5}$), and the magnitude of the
capital city, Budapest ($10^{6}$). In 2003, there is a minimum in the curves
of all size regimes, meaning that mid-sized source settlements act less as
invitors in the network than small villages and large towns. The contribution
to the $\log$ invitation count is highest for large towns and the capital, and
the relationship is significantly nonlinear between $\log w_{ST}$ and $\log
P_{s}$ (as confirmed by the coefficients of second-order terms in SI3). This
suggests that in the beginning of the life cycle of a product, the effect of
source settlement population size is even larger than what an urban scaling
law would predict. In 2006, $\log w_{ST}$ increases with $\log P_{S}$ in all
four selected $P_{T}$ regimes, therefore, the larger the source settlement,
the more the invitation count for a given target size. In 2011, this
monotonicity still holds, although the relationship becomes concave suggesting
a faster-than-power-law decrease in the role of large settlements in the
diffusion of invitations. For given source settlement sizes (again $10^{3}$,
$10^{4}$, $10^{5}$, and $10^{6}$), panels d-f show that an increasing target
settlement size increases invitation counts. The curves from a-f (as well as
coefficients in SI4) show that second-order terms are indeed important in
estimating the number of invitations, therefore, the observed relationships
are rarely linear between $\log w_{ST}$ and $\log P_{S}$ or $\log P_{T}$. This
means that power-laws between $\log w_{ST}$ and any combination of $P_{S}$ and
$P_{T}$ are not feasible models for this system.

But Figure 2 does not tell us about the effect of hierarchy, since we cannot
compare the combined effect and relative position of $P_{S}$ and $P_{T}$ on
$\log w_{ST}$. Panels a-c on Figure 3 present the level curves of the
contribution to $\log w_{ST}$ from the population terms and the constant $C$
from equation (1) in the variable space of $P_{S}$ and $P_{T}$. Coloring of
the level curves goes from the darkest minimum to the lightest maximum. In the
bottom left corner of panels a-c, $P_{S}$ and $P_{T}$ both correspond to small
village population sizes, whereas the top right corner corresponds to the case
where both $P_{S}$ and $P_{T}$ are as large as Budapest, the capital city. In
between, this area represents every possible (even hypothetical) settlement
pairs in Hungary in terms of $P_{S}$ and $P_{T}$. As expected, the number of
invitations is highest in all three selected years in the top right corner: if
both the source and the target settlements are very big. But the finer
structure of this estimation is remarkable, and is able to tell us more about
the hierarchy effects in detail.

Figure 3: Urban hierarchy in spatial diffusion. Contour lines in panels (a-c)
show level curves of the contribution of the source ($P_{S}$) and target
($P_{T}$) settlement sizes to the negative binomial term of the regression
with the distance $d$ fixed. Lighter color indicates more invitations, see the
colorbar of the top left corner. Green dashed lines correspond to
$(P_{S},P_{T})$ pairs where the number of potential total connections between
two settlements is constant, that is, where $P_{S}\cdot P_{T}$ is constant,
with darker color referring to more potential total connections ($10^{4}$,
$10^{6}$, $10^{8}$, or $10^{1}0$). Panels (d-f) show the surface height
corresponding to the position of the green dashed lines from panels (a-c).
These are the estimated invitation count contributions given that the
settlement pair has the same number of potential total contacts and fixed
distance.

As a most important marker of hierarchy effects, we investigate four different
cross-sections of this surface estimated by ZINB that are marked by the green
dashed lines of tangent -1 in the upper panels. These green dashed lines
correspond to $P_{S}\cdot P_{T}=\mathrm{const.}$ pairs, in other words, these
lines mark all settlement pairs for which the potential number of all
connections is the same. Along this line in panels a-c, in the very middle,
there is the case when $P_{S}=P_{T}$, but towards the right bottom corner,
$P_{S}>P_{T}$, and the population share of the source settlement grows.
Towards the top left corner, $P_{T}>P_{S}$, and the population share of the
target settlement grows for the fixed number of total connections. The bottom
panels d-f show the surface height along these selected lines, whose shape can
be classified into four different categories.

  1. 1.

If it is increasing towards large $P_{S}/P_{T}$ values, it means that the
invitation contribution is larger for those pairs, in which the source
settlement is larger than the target, even if there are the same number of
total possible connections for other pairs. This corresponds to the
hierarchical invitation pattern in which larger settlements tend to invite
smaller settlements.

  2. 2.

On the other hand, if this curve is monotonically decreasing, then it is
exactly the other way round: relatively larger targets and smaller sources
lead to more invitations. This corresponds to a reverse effect in which
invitations flow back from smaller settlements to larger ones.

  3. 3.

If there is a minimum in this curve within the given realistic size regimes,
it means that larger size differences between $P_{S}/P_{T}$ are favored in
both directions compared to settlements of roughly the same sizes, or that
large size differences in either direction lead to more invitations.

  4. 4.

A maximum would select one “favored” size difference for which invitation
count contribution is maximal.

We can see that the largest size regime always falls into the first category
of the above four, here, a relatively larger source settlement in the pair
yields more invitations (Case 1). However, in 2003, the two smallest
categories, $P_{S}\cdot P_{T}=10^{6}$, and $P_{S}\cdot P_{T}=10^{8}$ both show
a strong inclination for reverse hierarchical invitations (Case 2), where the
smaller $P_{S}/P_{T}$ is, the larger the invitation count. And the regime of
cities with roughly 100,000 population disfavors similar sizes: contribution
is largest at both ends of the possible $P_{S}/P_{T}$ regime (Case 4). In
2006, these remarkably diverging behaviours come close to each other, and in
the two largest size regimes, there is a small hierarchical effect (Case 1),
the third line has a very shallow minimum, meaning that hierarchy plays almost
no role here, and there is a very slight decreasing effect for the smallest
settlement pairs (Case 2). In 2011, the three largest curves show a favor for
similar settlement sizes (Case 3), and a slight hierarchical effect for the
smaller settlement pairs (Case 1).

Figure 4: Distance effect in spatial diffusion. Value of the coefficient
$\chi$ estimated from the negative binomial part of the ZINB model
characterizing the dependence of invitations on spatial distance between two
settlements.

Finally, Figure 4 shows the coefficient $\chi$ of the distance term $\log
d_{ST}$ from equation (1) with its estimated error. There is a clear
distinction between years with more long-range invitations (2002-2005,
2011-2012), and the short-range invitation years (2006-2010). This result from
the regression is thus in line with the intuitive picture from Figure 1, in
which the share of invitations sent to less than 20 km was very high in Figure
1c for most the most probable invitation tree edges in years 2006-2009. Here
we can see that this effect remains if we control for excessive zero counts
and population sizes.

##  3 Discussion

This paper contributes to the traditional discussion on spatial diffusion of
innovations by disentangling the role of urban hierarchy from geographical
distance using data from unprecedented detail and timescale. Looking at the
most probable diffusion paths using direction and timestamp of nearly 3
million spreading events, we discover that diffusion initially cascades from
large cities to all other settlements. Later, in the middle of the life cycle
when adoption accelerates, hierarchical diffusion manifests in long chains and
emerges along the entire spectrum of urban hierarchy. Unlike in previous
research that investigated distance effect separately, our multivariate
approach unfolds the real effect of geographical distance that intensifies
stronger around the adoption peak but becomes weak again in the late phases of
the life cycle.

Our results shed new light on the main puzzles of spatial diffusion research:
how is the spatial distribution of individuals open to novelty and how does
the structure of social interaction networks influence diffusion prediction at
local scales? Previous models have found that settlement size and distance are
remaining sources of error in complex contagion models. The finding that the
role of urban hierarchy is dynamically changing over the life cycle and has a
consequence on the role of distance as well might be important for future
models. Simulated complex contagion models on social networks should
incorporate the hierarchical diffusion processes.

The consequences of these findings might reach far concerning the potentials
of policy interventions in supporting or blocking diffusion. Adoption of new
technologies or scientific knowledge is the engine of progress; thus, a policy
targeting adoption according to the place of the settlement in the urban
hierarchy can benefit the society. However, the costs and motivations of
adoption vary greatly across products, and technologies; hence, the generality
of spatial patterns reported here must be looked at across various innovation
flows. The diffusion of social media is a specific example in which adoption
is greatly motivated by social interactions on the website. To prove the
generality of the illustrated spatial diffusion patterns, one needs to
investigate products that are not strongly associated with the social
dimension. Policy might want to block the diffusion of harmful novelty, such
as fake news. However, these typically spread very quickly; thus, future
research shall compare spatial diffusion of innovations of long and short life
cycles. Finally, other socio-economic factors such as development,
instrastructure, education, institutions among others might have a role in
diffusion that we did not consider in this paper.

##  4 Materials and methods

###  4.1 Data

Although free registration became possible in 2012, this previously singular
option remained the major means of joining the website, with more than 60% of
iWiW users registering through an invitor even after this date. Thus, for the
majority of the users, we have the anonymized identifier of the inviting user
as well. Since our data was collected in 2013, the ID of the invitor was
missing in those cases when the profile of invitation sender was already
deleted or the profile was registered after June 2012 via the free
registration option.

Users had to choose a town of residence to be displayed along their profiles.
Users could change this location field easily, moreover, the reliability of
this piece of information cannot be checked, although the menu offered a
limited number of valid location strings for users. Nevertheless, we use this
town of residence as geolocation information about both the users and their
invitors.

In our network, nodes are geolocated at the settlement level using the self-
reported location field, we leave users having more that one, but less than
5000 connections, and we omit users for whom the first and last login happened
on the very same day. Moreover, we only include users with a Hungarian
location, since a certain number of the registrations came from neighboring or
other European countries.

Individual data access is restricted by an NDA between the research group and
the data owner. Aggregated data can be requested in email sent to
<EMAIL_ADDRESS>
###  4.2 Most probable invitation paths

We construct the most probable spreading of the invitations on the settlement
networks by looking for the most probable invitation tree using an efficient
implementation of the Chu-Liu-Edmonds algorithm [31]. We suppose that the
probability of an invitation path is the product of the probabilities along
the constituting directed edges, where we normalize invitation counts by the
number of total invitations going out of the source settlement in three time
periods (2003-2005, 2006-2009, 2010-2012). By taking the logarithm and
multipying edge weights by -1, this problem is equivalent to finding the
minimum weight branching in the directed weighted invitation network. Because
the network is not always strongly connected, we add an artificial source node
that we connect to every other node with a sufficiently large edge weight such
that the algorithm avoids these edges if possible, but we remove the
artificial root node at the end of the process.

###  4.3 Regression framework

The zero-inflated negative binomial regression model is formalized by:

$\displaystyle w_{ST}$
$\displaystyle=e^{k}\cdot\underbrace{\left(1-\frac{e^{l}}{1+e^{l}}\right)}_{\mathrm{nonzero~{}inv.~{}probability}},$
(2) $\displaystyle\log w_{ST}$
$\displaystyle=k+\underbrace{\log\left(1-\frac{e^{l}}{1+e^{l}}\right)}_{\mathrm{small}}.$
(3)

In the first of the above equations (2), $l$ is the log-odds term from the
logit model that is used to estimate the propability that there is a nonzero
number of invitations between settlement pairs. The parameter $c$
characterizes the actual invitation counts between settlement pairs given that
the count is nonzero. The second equation (3) shows that the contribution to
$\log w_{ST}$ is mostly contained in the term $k$, since the probability of
the invitation between the two settlements being nonzero is close to one, for
which the logarithm is close to 0. Both $k$ and $l$ are governed by the
processes mentioned before, therefore, we choose the terms that predict them
in the ZINB model to reflect these.

We estimate the terms $k$ and $l$ from Equation 2 for each year between
2002-2012 using Stata [34].

##  5 Acknowledgements

The suggestions of Ferenc Bencs in identifying most probable invitation paths
is gratefully acknowledged. Figure 1 has been finalized by graphic designer
Szabolcs Tóth-Zs. The authors acknowledge financial help received from
National Research, Development and Innovation Office of Hungary (KH 130502).

## References

  * [1] Joseph Schumpeter and Ursula Backhaus  “The theory of economic development”  In _Joseph Alois Schumpeter_ Springer, 2003, pp. 61–116
  * [2] Robert M Solow  “A contribution to the theory of economic growth”  In _Quarterly Journal of Economics_ 70.1 MIT Press, 1956, pp. 65–94  DOI: 10.2307/1884513
  * [3] Frank M Bass  “A new product growth for model consumer durables”  In _Management science_ 15.5 INFORMS, 1969, pp. 215–227
  * [4] Everett M Rogers  “Diffusion of innovations”  SimonSchuster, 2010
  * [5] Zvi Griliches  “Hybrid Corn: An Exploration in the Economics of Technological Change”  In _Econometrica_ 25.4 JSTOR, 1957, pp. 501  DOI: 10.2307/1905380
  * [6] Torsten Hägerstrand  “Innovation Diffusion as a Spatial Process”  Chicago, USA: Univ. Chicago Press., 1967
  * [7] Peter R Gould  “Spatial diffusion: the spread of ideas and innovations in geographic space”  Learning Resources in International Studies, 1975
  * [8] Ping Luo et al.  “A study on cellular automata based on relational databases and spatio-temporal simulations of culture diffusion”  In _Chinese Geographical Science_ 12.4, 2002, pp. 359–365  DOI: 10.1007/s11769-002-0043-7
  * [9] S Bauernschuster and O Falck  “Culture, spatial diffusion of ideas and their long-lasting imprints–evidence from Froebel’s kindergarten movement”  In _Journal of Economic Geography_ 15.3, 2015, pp. 601–630  DOI: 10.1093/jeg/lbu028
  * [10] Sendy Farag et al.  “E-Shopping in the Netherlands: Does Geography Matter?”  In _Environment and Planning B: Planning and Design_ 33.1, 2006, pp. 59–74  DOI: 10.1068/b31083
  * [11] Jure Leskovec and Eric Horvitz  “Geospatial structure of a planetary-scale social network”  In _IEEE Transactions on Computational Social Systems_ 1.3 IEEE, 2014, pp. 156–163  DOI: 10.1109/TCSS.2014.2377789
  * [12] Walter Christaller  “Central places in southern Germany”  Prentice-Hall, 1966
  * [13] David Strang and Nancy Brandon Tuma  “Spatial and Temporal Heterogeneity in Diffusion”  In _American Journal of Sociology_ 99.3 University of Chicago Press, 1993, pp. 614–639  DOI: 10.1086/230318
  * [14] Thomas W Valente  “Social network thresholds in the diffusion of innovations”  In _Social Networks_ 18.1 Elsevier, 1996, pp. 69–89  DOI: 10.1016/0378-8733(95)00256-1
  * [15] Damon Centola and Michael Macy  “Complex contagions and the weakness of long ties”  In _American Journal of Sociology_ 113.3 The University of Chicago Press, 2007, pp. 702–734  DOI: 10.1086/521848
  * [16] Damon Centola  “The spread of behavior in an online social network experiment”  In _Science_ 329.5996 American Association for the Advancement of Science, 2010, pp. 1194–1197  DOI: 10.1126/science.1185231
  * [17] Marton Karsai et al.  “Local cascades induced global contagion: How heterogeneous thresholds, exogenous effects, and unconcerned behaviour govern online adoption spreading”  In _Scientific Reports_ 6 Nature Publishing Group, 2016  DOI: 10.1038/srep27178
  * [18] Johan Ugander, Lars Backstrom, Cameron Marlow and Jon Kleinberg  “Structural diversity in social contagion”  In _Proceedings of the National Academy of Sciences of the United States of America_ 109.16 National Acad Sciences, 2012, pp. 5962–5966  DOI: 10.1073/pnas.1116502109
  * [19] Jameson L Toole, Meeyoung Cha and Marta C González  “Modeling the adoption of innovations in the presence of geographic and media influences”  In _PLoS ONE_ 7.1 Public Library of Science, 2012, pp. e29528  DOI: 10.1371/journal.pone.0029528
  * [20] Balázs Lengyel et al.  “The role of geography in the complex diffusion of innovations”  In _Scientific Reports_ 10.15065 Nature Research, 2020  DOI: 10.1038/s41598-020-72137-w
  * [21] Michael Trusov, William Rand and Yogesh V Joshi  “Improving prelaunch diffusion forecasts: Using synthetic networks as simulated priors”  In _Journal of Marketing Research_ 50.6 American Marketing Association, 2013, pp. 675–690  DOI: 10.1509/jmr.11.0508
  * [22] Balázs Lengyel et al.  “Geographies of an online social network”  In _PLoS ONE_ 10.9 Public Library of Science, 2015, pp. e0137248  DOI: 10.1371/journal.pone.0137248
  * [23] Johannes Wachs, Taha Yasseri, Balázs Lengyel and János Kertész  “Social capital predicts corruption risk in towns”  In _Royal Society Open Science_ 6.4 The Royal Society, 2019, pp. 182103  DOI: 10.1098/rsos.182103
  * [24] Gergő Tóth et al.  “Inequality is rising where social network segregation interacts with urban topology”  In _Nature Communications_ 12.1, 2021, pp. 1143  DOI: 10.1038/s41467-021-21465-0
  * [25] János Török and János Kertész  “Cascading collapse of online social networks”  In _Scientific Reports_ 7.1 Nature Publishing Group, 2017, pp. 16743  DOI: 10.1038/s41598-017-17135-1
  * [26] László Lőrincz, Júlia Koltai, Anna Fruzsina Győr and Károly Takács  “Collapse of an online social network: Burning social capital to create it?”  In _Social Networks_ 57 Elsevier, 2019, pp. 43–53  DOI: 10.1016/j.socnet.2018.11.004
  * [27] Zsolt Katona, Peter Pal Zubcsek and Miklos Sarvary  “Network effects and personal influences: The diffusion of an online social network”  In _Journal of Marketing Research_ 48.3 American Marketing Association, 2011, pp. 425–443  DOI: 10.1509/jmkr.48.3.425
  * [28] Ana Pastore y Piontti et al.  “The infection tree of global epidemics”  In _Network Science_ 2.1, 2014, pp. 132–137  DOI: 10.1017/nws.2014.5
  * [29] Frank Schlosser et al.  “COVID-19 lockdown induces disease-mitigating structural changes in mobility networks”  In _Proceedings of the National Academy of Sciences of the United States of America_ 117.52 National Acad Sciences, 2021, pp. 32883–32890  DOI: 10.1073/PNAS.2012326117
  * [30] Robert Endre Tarjan  “Finding optimum branchings”  In _Networks_ 7.1 Wiley Online Library, 1977, pp. 25–35  DOI: 10.1002/net.3230070103
  * [31] Harold N Gabow, Zvi Galil, Thomas Spencer and Robert E Tarjan  “Efficient algorithms for finding minimum spanning trees in undirected and directed graphs”  In _Combinatorica_ 6.2, 1986, pp. 109–122  DOI: 10.1007/BF02579168
  * [32] Luís M.A. Bettencourt et al.  “Growth, innovation, scaling, and the pace of life in cities”  In _Proceedings of the National Academy of Sciences of the United States of America_ 104.17 National Acad Sciences, 2007, pp. 7301–7306  DOI: 10.1073/pnas.0610172104
  * [33] Scott J. Long  “Regression Models for Categorical and Limited Dependent Variables”, 1997, pp. 328  URL: https://us-sagepub-com.vu-nl.idm.oclc.org/en-us/nam/regression-models-for-categorical-and-limited-dependent-variables/book6071
  * [34] “Stata 17 help for zinb”  URL: https://www.stata.com/help.cgi?zinb

## Supplementary Information

## Supplementary Information 1: The iWiW life cycle

The number of registrations was low in the first few years but iWiW reached a
several hundred thousands of people by 2005. The largest Hungarian internet
provider company acquired the website in 2006, and the number of registered
users grew rapidly from 1.5 to more than 4 millions until December 2008.
Moreover, between 2005 and 2010, iWiW was the most frequently visited
Hungarian webpage. This is well reflected in Figure 5, that shows the number
of registrations over the full life cycle of the website, and confirms that at
the peak of iWiW’s popularity, the number of registered users jumped to more
than 50.000 monthly, and increased further with high variability to a peak
about 90.000 invitations per month until the middle of 2007.

Figure 5: Number of registrations on the iWiW social network by month. Red
curve shows the cumulative number of adopters, annotations reflect the adopter
categories of [4].

The figure shows how early adoption was driven by the age group of 19-29
years, then people of age 8-18 and 30-40 followed. Older age groups between
41-51, 52-62, then 63-73 gradually joined in the adoption process. This means
that age structure and age-related social and family contacts of settlements
might influence adoption patterns.

Figure 6: Number of registrations on the iWiW social network by month per age
group.

## Supplementary Information 2: The share of invitations by diffusion
categories

The average distance of invitations declines somewhat in 2006-2007, although
the standard deviation of invitation distances marked by the errorbars remains
high. However, this distance decrease is not fully in line with the regression
coefficient of the $\log d$ term, because here, we do not control for
hierarchical patterns.

Figure 7: Average distance of invitations in time.

The average population size difference between source and target settlement
decreases, that might be connected to the decreasing role of hierarchy in the
invitation process.

Figure 8: Average population difference of invitations in time.

High share of short-range invitations in the middle of the life cycle is in
line with the regression coefficient of the distance term from the negative
binomial part of the ZINB regression.

Figure 9: Share of invitations sent to less than 20 km.

The share of invitations sent to a settlement at least 3x smaller than the
source settlement is constantly decreasing. In the early and late phases,
share of invitations sent into the different direction, to an at least 3x
larger settlement than the source settlement is somewhat larger than in the
middle of the time range.

Figure 10: Share of invitations downwards and upwards the hierarchy level (at
least 3x population difference in either direction).

The yearly distance distributions show a very pronounced tail for the early
years indicating the large weight of long-range invitations in the early
stages. High share of short-range invitations indicate the presence of
neighborhood diffusion.

Figure 11: Share of invitations in given distance ranges per year.

The distribution of hierarchy differences in invitations shows a tendency
towards hierarchical diffusion, since there are always more invitations sent
to smaller settlements (x axis larger than 0), than invitation sent from
smaller to larger settlements (x axis smaller than 0). In the early stages,
most invitations go between settlements of very different sizes, then similar
hierarchy levels tend to dominate the diffusion process.

Figure 12: Share of invitations with given $log_{10}$ population difference
ranges per year.

## Supplementary Information 3: Regression coefficients

In the following two tables, we can find the regression coefficients of the
full zero-inflated negative binomial model for all years. Significance levels
are marked by stars, standard errors of coefficients are given in brackets.
The inflation part corresponds to the estimation of the logit probability term
$l$ from equation (2), and the negative binomial part to the estimation of the
count process term $k$ from (2). Both terms are estimated by the form

$\alpha(\log P_{S})^{2}+\beta\cdot(\log P_{S}\cdot\log P_{T})+\gamma(\log
P_{T})^{2}+\delta\log P_{S}+\varepsilon\log P_{T}+\chi\log d_{ST}+C,$

where $\alpha$, $\beta$, $\gamma$, $\delta$, $\varepsilon$, and $\chi$ are the
coefficients of the variables listed in the regression tables. Note the
significance of the second-order terms in the negative binomial part of the
model.

| Year | 2002 | 2003 | 2004 | 2005 | 2006 | 2007
---|---|---|---|---|---|---|---
|  | Coeff. | SE | Coeff. | SE | Coeff. | SE | Coeff. | SE | Coeff. | SE | Coeff. | SE
Model part | Variable |  |  |  |  |  |  |  |  |  |  |  |
| Num. obs. | 6525466 |  | 6525466 |  | 6525466 |  | 6525466 |  | 6525466 |  | 6525466 |
| $\log\alpha$ | -0.568*** | (0.176) | -0.832*** | (0.147) | -0.895*** | (0.126) | -0.792*** | (0.0511) | -0.274*** | (0.0156) | -0.0661*** | (0.0121)
Inflation | Constant | 28.52*** | (10.84) | 32.24*** | (5.794) | 23.99*** | (3.397) | 13.05*** | (1.433) | -7.983*** | (0.78) | -10.57*** | (0.665)
| $\log P_{S}$ | -2.452 | (1.514) | -2.531*** | (0.629) | -2.012*** | (0.339) | -0.65*** | (0.178) | 1.298*** | (0.103) | 1.635*** | (0.0974)
| $\log P_{S}^{2}$ | 0.00491 | (0.102) | 0.0241 | (0.0254) | 0.023* | (0.0119) | -0.026*** | (0.00704) | -0.111*** | (0.00457) | -0.129*** | (0.00512)
| $\log P_{S}\cdot\log P_{T}$ | 0.0686 | (0.0473) | 0.0709*** | (0.023) | 0.039*** | (0.0151) | 0.00279 | (0.00816) | -0.0171*** | (0.00622) | -0.0309*** | (0.00636)
| $\log P_{T}$ | -1.777** | (0.845) | -2.201*** | (0.55) | -1.215*** | (0.369) | -0.69*** | (0.157) | 1.515*** | (0.107) | 1.62*** | (0.0945)
| $\log P_{T}^{2}$ | 0.0417 | (0.0296) | 0.0434*** | (0.0165) | 0.00774 | (0.0115) | -0.00393 | (0.00505) | -0.119*** | (0.00513) | -0.114*** | (0.0045)
| $\log d_{ST}$ | 0.383*** | (0.122) | 0.487*** | (0.071) | 0.625*** | (0.0511) | 0.68*** | (0.0236) | 0.6*** | (0.0124) | 0.63*** | (0.0105)
Neg. binom. | Constant | 15.55 | (10.94) | 22.22*** | (5.611) | 17.81*** | (3.306) | 9.638*** | (1.201) | -1.285*** | (0.458) | -2.62*** | (0.295)
| $\log P_{S}$ | -4.256*** | (0.982) | -4.205*** | (0.47) | -3.575*** | (0.301) | -2.178*** | (0.122) | -0.532*** | (0.0514) | -0.29*** | (0.0355)
| $\log P_{S}^{2}$ | 0.21*** | (0.0231) | 0.176*** | (0.0114) | 0.145*** | (0.00856) | 0.0955*** | (0.00358) | 0.0506*** | (0.0016) | 0.0459*** | (0.00122)
| $\log P_{S}\cdot\log P_{T}$ | 0.06 | (0.0414) | 0.115*** | (0.0218) | 0.113*** | (0.0128) | 0.104*** | (0.00536) | 0.0625*** | (0.00252) | 0.0416*** | (0.00186)
| $\log P_{T}$ | -0.778 | (0.897) | -1.602*** | (0.519) | -1.143*** | (0.332) | -0.882*** | (0.125) | -0.0645 | (0.0511) | 0.31*** | (0.0354)
| $\log P_{T}^{2}$ | 0.0649*** | (0.0187) | 0.0661*** | (0.0124) | 0.0416*** | (0.00933) | 0.0369*** | (0.00358) | 0.0231*** | (0.00158) | 0.0123*** | (0.0012)
| $\log d_{ST}$ | -1.051*** | (0.0526) | -0.882*** | (0.039) | -0.787*** | (0.0342) | -0.918*** | (0.0161) | -1.389*** | (0.00789) | -1.605*** | (0.00659)
(* = $p<0.10$, **=$p<0.05$, ***=$p<0.01$)
Table 1: Results of the zero-inflated negative binomial regression models for
each year (2002-2007).

| Year | 2008 | 2009 | 2010 | 2011 | 2012 |  |
---|---|---|---|---|---|---|---|---
|  | Coeff. | SE | Coeff. | SE | Coeff. | SE | Coeff. | SE | Coeff. | SE |  |
Model part | Variable |  |  |  |  |  |  |  |  |  |  |  |
| Num. obs. | 6525466 |  | 6525466 |  | 6525466 |  | 6525466 |  | 6525466 |  |  |
| $\log\alpha$ | -0.123*** | (0.0171) | -0.249*** | (0.0218) | -0.315*** | (0.0248) | -1.24*** | (0.0788) | -2.817* | (1.494) |  |
Inflation | Constant | -7.207*** | (0.619) | -6.497*** | (0.644) | -6.925*** | (0.671) | 0.316 | (0.952) | 13.23*** | (3.106) |  |
| $\log P_{S}$ | 1.212*** | (0.0869) | 1.237*** | (0.0902) | 1.153*** | (0.0913) | 0.379*** | (0.122) | -1.171*** | (0.375) |  |
| $\log P_{S}^{2}$ | -0.094*** | (0.00409) | -0.106*** | (0.00424) | -0.1*** | (0.00419) | -0.0939*** | (0.0054) | -0.0482*** | (0.0173) |  |
| $\log P_{S}\cdot\log P_{T}$ | -0.0207*** | (0.00577) | 0.0077 | (0.00578) | 0.0124** | (0.00565) | 0.0812*** | (0.00666) | 0.132*** | (0.0189) |  |
| $\log P_{T}$ | 1.006*** | (0.0878) | 0.791*** | (0.0891) | 0.893*** | (0.0901) | 0.21* | (0.12) | -0.439 | (0.388) |  |
| $\log P_{T}^{2}$ | -0.0809*** | (0.00413) | -0.0821*** | (0.00418) | -0.0864*** | (0.00413) | -0.0858*** | (0.00531) | -0.0886*** | (0.0174) |  |
| $\log d_{ST}$ | 0.657*** | (0.0123) | 0.723*** | (0.0139) | 0.739*** | (0.0155) | 1.087*** | (0.0249) | 1.208*** | (0.0791) |  |
Neg. binom. | Constant | -1.812*** | (0.31) | -1.763*** | (0.349) | -2.677*** | (0.369) | -1.183 | (0.725) | 9.663*** | (3.017) |  |
| $\log P_{S}$ | 0.00134 | (0.0374) | 0.0804* | (0.0416) | 0.213*** | (0.0438) | 0.0565 | (0.0809) | -1.367*** | (0.318) |  |
| $\log P_{S}^{2}$ | 0.031*** | (0.00134) | 0.0222*** | (0.0015) | 0.0145*** | (0.00158) | -0.00572** | (0.00285) | 0.0211* | (0.0116) |  |
| $\log P_{S}\cdot\log P_{T}$ | 0.034*** | (0.00194) | 0.0386*** | (0.0021) | 0.0362*** | (0.00219) | 0.0748*** | (0.00358) | 0.141*** | (0.0114) |  |
| $\log P_{T}$ | -0.096** | (0.0387) | -0.223*** | (0.0428) | -0.144*** | (0.0445) | -0.442*** | (0.0805) | -1.447*** | (0.322) |  |
| $\log P_{T}^{2}$ | 0.0301*** | (0.00139) | 0.0318*** | (0.00154) | 0.0275*** | (0.0016) | 0.0143*** | (0.00279) | 0.0218* | (0.0117) |  |
| $\log d_{ST}$ | -1.566*** | (0.00909) | -1.476*** | (0.0111) | -1.442*** | (0.0131) | -1.017*** | (0.0212) | -0.709*** | (0.066) |  |
(* = $p<0.10$, **=$p<0.05$, ***=$p<0.01$)
Table 2: Results of the zero-inflated negative binomial regression models for
each year (2007-2012).