about $\sigma^{d-1}$. Here, choosing the radius suitably means taking the
radius to be the closest value in $E$ to the $x_{d}$ coordinate and then
rescaling by $\frac{1}{\sqrt{2}}$.
We notice that $R$ has volume about $(\delta/\sigma)^{d-1}\delta N(E\cap
I,\delta)$. Putting the lower bound (89) into (80) leads to the inequality

$\sigma^{2d-2}\delta\left((\delta/\sigma)^{d-1}\delta N(E\cap
I,\delta)\right)^{1/r}\leq(\sigma^{d-1}\delta)^{p^{-1}+q^{-1}}$ (90)

Recalling that $\sigma=\delta^{\alpha/2}$, the claim follows by comparing the
exponents and doing routine algebraic manipulation. ∎

Putting Propositions 7.1, 7.3, and 7.5 all together we get the following
necessary conditions for the boundedness of $\mathcal{A}_{E}$.

######  Proposition 7.6 (Necessary conditions for boundedness of
$\mathcal{A}_{E}$).

Let $E\subset[1,2]$ be a Assouad regular set with $\text{dim}_{M}(E)=\beta$
and $\text{dim}_{A}(E)=\gamma$. If $\mathcal{A}_{E}:L^{p}\times
L^{q}\rightarrow L^{r}$ is bounded, then

$\frac{1}{r}\leq\frac{1}{p}+\frac{1}{q}\leq 1+m_{linear}(d,r,\beta,\gamma).$

where

$m_{linear}(d,r,\beta,\gamma):=\text{min}\left\\{\frac{d-1}{d}+\frac{1-\beta}{dr},\frac{\beta(d-1)}{\beta(d-1)+2\gamma}+\frac{(d-\beta)2\gamma-(d-1)\beta}{\beta(d-1)+2\gamma}\frac{1}{r},\frac{d}{r}\right\\}.$

The example of Proposition 7.1 can also be adapted to give necessary
conditions in the biparameter setting.

######  Proposition 7.7.

Suppose $E_{1},E_{2}$ have upper Minkowski dimensions $\beta_{1},\beta_{2}$.
Define the quantity

$\beta^{*}=\dim_{M}(E_{1}E_{2}(E_{1}^{2}+E_{2}^{2})^{-1/2}).$

If $\mathcal{A}^{(2)}_{S^{2d-1},E_{1},E_{2}}:L^{p}\times L^{q}\rightarrow
L^{r}$ continuously then we have

$\frac{1}{p}+\frac{1}{q}\leq\frac{2d-1}{d}+\frac{1-\beta^{*}}{dr}$ (91)

###### Proof.

Let $f=g=\chi_{B(0,C\delta)}$ and denote
$h=\mathcal{A}^{(2)}_{\hat{\sigma},E_{1},E_{2}}(f,g)$. Mimicking the single
parameter case, we want to find the set of $x$ where we have the lower bound
$|h(x)|\gtrsim\delta^{2d-1}$. Since we now have parameters $t_{1},t_{2}$, the
region where this lower bound holds contains all $x$ such that there exist
$t_{1}\in E_{1},t_{2}\in E_{2}$ and some $|y|\in(0,1)$ satisfying

$|x|=t_{1}|y|=t_{2}\sqrt{1-|y|^{2}}.$

Really, the first equation just needs to hold up to $\delta$ error so that we
can compute the $\delta$-covering number and take a limit. Substituting the
second equation necessitates that $|y|$ solves the quadratic equation

$\left(\frac{t_{1}^{2}+t_{2}^{2}}{t_{2}^{2}}\right)|y|^{2}=1$

and so we can write the solution as

$|x|=t_{1}|y|=\frac{t_{1}t_{2}}{\sqrt{t_{1}^{2}+t_{2}^{2}}}$

Thus, the values of $|x|$ that are admissible is a $C\delta$-neighborhood of
points of the form

$E^{*}:=\left\\{\frac{t_{1}t_{2}}{\sqrt{t_{1}^{2}+t_{2}^{2}}}:t_{1}\in
E_{1},t_{2}\in E_{2}\right\\}$

Call $\beta^{*}=\dim_{M}(E^{*})$, which depends on $E_{1},E_{2}$. Testing this
against the boundedness condition gives

$\delta^{2d-1}\delta^{1/r}N(E^{*},\delta)^{1/r}\leq\delta^{d/p}\delta^{d/q}$
(92)

which leads to the inequality

$\frac{d}{p}+\frac{d}{q}\leq 2d-1+\frac{1-\beta^{*}+\epsilon}{r}$

We can send $\epsilon$ to 0 in the limit, so dividing through by $d$ gives the
claim. ∎

Similarly, we can adapt Proposition 7.4 to the biparameter setting.

######  Proposition 7.8.

Suppose $E_{1},E_{2}$ have upper Minkowski dimension $\beta_{1},\beta_{2}$ and
let $\beta^{*}$ be as in the previous proposition. If
$\mathcal{A}_{E}:L^{p}\times L^{q}\rightarrow L^{r}$ continuously then we have

$\frac{1}{p}+\frac{1}{q}\leq\frac{2d}{d+1}+\frac{1}{r}\left(1-\frac{2\beta^{*}}{(d+1)}\right)$
(93)

###### Proof.

Once again, we can define

$R_{1}=[-C_{1}\sqrt{\delta},C_{1}\sqrt{\delta}]^{d-1}\times[-C_{1}\delta,C_{1}\delta],\quad
R_{2}=[-C_{2}\sqrt{\delta},C_{2}\sqrt{\delta}]^{d-1}\times[-C_{2}\delta,C_{2}\delta],$
$R_{3}=[-\sqrt{\delta},\sqrt{\delta}]^{d-1}\times\left(\bigcup_{i=1}^{N(E^{*},\delta)}I_{i}\right)$

where we may take $C_{1},C_{2}=100$, $\delta$ is sufficiently small, and
$\\{I_{i}\\}$ is a minimal collection of $\delta$-intervals covering the set
$E^{*}$. Then for $x\in R_{3}$, one can check that

$\mathcal{A}_{E}(f_{\delta},g_{\delta})(x)\gtrsim\delta^{d}$

That combined with $|R_{3}|\sim\delta^{\frac{d-1}{2}}N(E^{*},\delta)\delta$ is
enough to finish the computation. ∎

######  Remark 7.9.

In fact, Proposition 7.5 also generalizes to the biparameter setting, with
$\beta,\gamma$ being replaced by $\beta^{*}$ and $\gamma^{*}$ respectively,
the Minkowski and Assouad dimensions of the set $E^{*}$. The details are very
similar to the other biparameter arguments we have given, so we leave them to
the interested reader.

In fact, we can say more about $\beta^{*}$ by making use of Lemma 2.1. The
dimension of $E^{*}$ is equal to the dimension of
$\frac{E_{1}^{2}E_{2}^{2}}{E_{1}^{2}+E_{2}^{2}}$ which can be rewritten as
$(E_{1}^{-2}+E_{2}^{-2})^{-1}$. Notice that if we freeze one of the sets
$E_{1}$ or $E_{2}$ and vary the $t_{i}$s in the other set, we immediately see
from Lemma 2.1 that $\beta^{*}\geq\max\\{\beta_{1},\beta_{2}\\}$ since
dimension is preserved under taking reciprocals, squaring, and translation by
a constant. This lower bound can be achieved, such as in the case that
$\beta_{1}=1$. On the other hand, given two sets $S_{1},S_{2}$, their sumset
satisfies the dimension bound

$\dim_{M}(S_{1}+S_{2})\leq\dim_{M}(S_{1})+\dim_{M}(S_{2}).$

Thus, the description

$\max\\{\beta_{1},\beta_{2}\\}\leq\beta^{*}\leq\beta_{1}+\beta_{2}$

is the best we can say for an arbitrary pair of sets
$E_{1},E_{2}\subset[1,2]$.

## References

  * [AHRS21] T. Anderson, K. Hughes, J. Roos, and A. Seeger.  $L^{p}\to L^{q}$ bounds for spherical maximal operators.  Math. Z., 297(3-4):1057–1074, 2021.
  * [BCAOR19] Alex Barron, José M. Conde-Alonso, Yumeng Ou, and Guillermo Rey.  Sparse domination and the strong maximal function.  Advances in Mathematics, 345:1–26, 2019.
  * [BCSS23] Ankit Bhojak, Surjeet Singh Choudhary, Saurabh Shrivastava, and Kalachand Shuin.  Sharp endpoint $L^{p}-$estimates for bilinear spherical maximal functions.  preprint arXiv:2310.00425, 2023.
  * [BF23] Tainara Borges and Benjamin Foster.  Bounds for lacunary bilinear spherical and triangle maximal functions.  preprint arXiv:2305.12269, 2023.
  * [BFO+23] Tainara Borges, Benjamin Foster, Yumeng Ou, Jill Pipher, and Zirui Zhou.  Sparse bounds for the bilinear spherical maximal function.  Journal of the London Mathematical Society, 107(4):1409–1449, 2023\.
  * [BGH+18] J. A. Barrionuevo, Loukas Grafakos, Danqing He, Petr Honzík, and Lucas Oliveira.  Bilinear spherical maximal function.  Math. Res. Lett., 25(5):1369–1388, 2018.
  * [Bou85] Jean Bourgain.  Estimations de certaines fonctions maximales.  C. R. Acad. Sci. Paris Sér. I Math., 301(10):499–502, 1985\.
  * [BRS20] David Beltran, Joris Roos, and Andreas Seeger.  Multi-scale sparse domination.  Mem. Amer. Math. Soc., to appear, preprint arXiv:2009.00022, 2020\.
  * [BS98] Jong-Guk Bak and Yong-Sun Shim.  Endpoint inequalities for spherical multilinear convolutions.  J. Funct. Anal., 157(2):534–553, 1998.
  * [CADPPV22] José M. Conde-Alonso, Francesco Di Plinio, Ioannis Parissis, and Manasa N. Vempati.  A metric approach to sparse domination.  Ann. Mat. Pura Appl. (4), 201(4):1639–1675, 2022.
  * [Cal79] Calixto P. Calderón.  Lacunary spherical means.  Illinois J. Math., 23(3):476–484, 1979.
  * [CDPO18] Amalia Culiuc, Francesco Di Plinio, and Yumeng Ou.  Domination of multilinear singular integrals by positive sparse forms.  J. Lond. Math. Soc. (2), 98(2):369–392, 2018.
  * [CGH+22] Jiecheng Chen, Loukas Grafakos, Danqing He, Petr Honzík, and Lenka Slavíková.  Bilinear maximal functions associated with surfaces.  Proc. Amer. Math. Soc., 150(4):1635–1639, 2022.
  * [CLS23] Chu-hee Cho, Jin Bong Lee, and Kalachand Shuin.  $L^{p}$ improving properties and maximal estimates for certain multilinear averaging operators.  preprint arXiv:2306.04196, 2023.
  * [CO18] Laura Cladek and Yumeng Ou.  Sparse domination of Hilbert transforms along curves.  Math. Res. Lett., 25(2):415–436, 2018.
  * [CW78] R. R. Coifman and Guido Weiss.  Review: R. E. Edwards and G. I. Gaudry, Littlewood-Paley and multiplier theory.  Bulletin of the American Mathematical Society, 84(2):242 – 250, 1978.
  * [CZ22] Michael Christ and Zirui Zhou.  A class of singular bilinear maximal functions, 2022.
  * [DR22] Georgios Dosidis and João P. G. Ramos.  The multilinear spherical maximal function in one dimension, 2022.
  * [DV98] Javier Duoandikoetxea and Ana Vargas.  Maximal operators associated to Fourier multipliers with an arbitrary set of parameters.  Proc. Roy. Soc. Edinburgh Sect. A, 128(4):683–696, 1998.
  * [GHH21] Loukas Grafakos, Danqing He, and Petr Honzík.  Maximal operators associated with bilinear multipliers of limited decay.  J. Anal. Math., 143(1):231–251, 2021.
  * [GHS20] Loukas Grafakos, Danqing He, and Lenka Slavíková.  $L^{2}\times L^{2}\to L^{1}$ boundedness criteria.  Math. Ann., 376(1-2):431–455, 2020.
  * [GPP23] José Gaitan, Eyvindur Ari Palsson, and Georgios Psaromiligkos.  On restricted falconer distance sets.  preprint arXiv:2305.18053, 2023.
  * [Gra14] Loukas Grafakos.  Modern Fourier analysis, volume 250 of Graduate Texts in Mathematics.  Springer, New York, third edition, 2014.
  * [HHY20] Yaryong Heo, Sunggeum Hong, and Chan Woo Yang.  Improved bounds for the bilinear spherical maximal operators.  Math. Res. Lett., 27(2):397–434, 2020.
  * [IPS22] Alex Iosevich, Eyvindur Ari Palsson, and Sean R. Sovine.  Simplex averaging operators: quasi-Banach and $L^{p}$-improving bounds in lower dimensions.  J. Geom. Anal., 32(3):Paper No. 87, 16, 2022.
  * [JL20] Eunhee Jeong and Sanghyuk Lee.  Maximal estimates for the bilinear spherical averages and the bilinear Bochner-Riesz operators.  J. Funct. Anal., 279(7):108629, 29, 2020.
  * [Lac17] Michael T. Lacey.  An elementary proof of the $A_{2}$ bound.  Israel J. Math., 217(1):181–195, 2017.
  * [Lac19] Michael T. Lacey.  Sparse bounds for spherical maximal functions.  J. Anal. Math., 139(2):613–635, 2019.
  * [Ler13] Andrei K. Lerner.  A simple proof of the $A_{2}$ conjecture.  Int. Math. Res. Not. IMRN, (14):3159–3170, 2013.
  * [Lit63] Walter Littman.  Fourier transforms of surface-carried measures and differentiability of surface averages.  Bull. Amer. Math. Soc., 69:766–770, 1963.
  * [Lit73] Walter Littman.  $L^{p}-L^{q}$-estimates for singular integral operators arising from hyperbolic equations.  In Partial differential equations (Proc. Sympos. Pure Math., Vol. XXIII, Univ. California, Berkeley, Calif., 1971), volume Vol. XXIII of Proc. Sympos. Pure Math., pages 479–481. Amer. Math. Soc., Providence, RI, 1973.
  * [Obe88] Daniel M. Oberlin.  Multilinear convolutions defined by measures on spheres.  Trans. Amer. Math. Soc., 310(2):821–835, 1988.
  * [PS22] Eyvindur Ari Palsson and Sean R. Sovine.  Sparse bounds for maximal triangle and bilinear spherical averaging operators.  preprint arXiv:2110.08928, 2022.
  * [PS23] Eyvindur A. Palsson and Sean R. Sovine.  The triangle operator.  preprint arXiv:1910.01282, 2023.
  * [RdF86] José L. Rubio de Francia.  Maximal functions and Fourier transforms.  Duke Math. J., 53(2):395–404, 1986.
  * [RS23] Joris Roos and Andreas Seeger.  Spherical maximal functions and fractal dimensions of dilation sets.  Amer. J. Math., 145(4):1077–1110, 2023.
  * [RSS21] Luz Roncal, Saurabh Shrivastava, and Kalachand Shuin.  Bilinear spherical maximal functions of product type.  J. Fourier Anal. Appl., 27(4):Paper No. 73, 42, 2021.
  * [Sch97] W. Schlag.  A generalization of Bourgain’s circular maximal theorem.  J. Amer. Math. Soc., 10(1):103–122, 1997.
  * [SS97] Wilhelm Schlag and Christopher D. Sogge.  Local smoothing estimates related to the circular maximal theorem.  Math. Res. Lett., 4(1):1–15, 1997.
  * [SS21] Saurabh Shrivastava and Kalachand Shuin.  Lp estimates for multilinear convolution operators defined with spherical measure.  Bulletin of the London Mathematical Society, 53(4):1045–1060, Mar 2021.
  * [Ste76] Elias M. Stein.  Maximal functions. I. Spherical means.  Proc. Nat. Acad. Sci. U.S.A., 73(7):2174–2175, 1976.
  * [STW03] Andreas Seeger, Terence Tao, and James Wright.  Endpoint mapping properties of spherical maximal operators.  Journal of the Institute of Mathematics of Jussieu, 2(1):109–144, 2003.
  * [SWW95] Andreas Seeger, Stephen Wainger, and James Wright.  Pointwise convergence of spherical means.  Math. Proc. Cambridge Philos. Soc., 118(1):115–124, 1995.
  * [SWW97] Andreas Seeger, Stephen Wainger, and James Wright.  Spherical maximal operators on radial functions.  Mathematische Nachrichten, 187(1):241–265, 1997.


# A Customizable Dynamic Scenario Modeling and Data Generation Platform for
Autonomous Driving

Jay Shenoy
UC Berkeley
<EMAIL_ADDRESS>
Edward Kim
UC Berkeley
<EMAIL_ADDRESS>
Xiangyu Yue
UC Berkeley
<EMAIL_ADDRESS>
Taesung Park
UC Berkeley
<EMAIL_ADDRESS>
Daniel Fremont
UC Santa Cruz
<EMAIL_ADDRESS>
Alberto Sangiovanni-Vincentelli
UC Berkeley
<EMAIL_ADDRESS>
Sanjit Seshia
UC Berkeley
<EMAIL_ADDRESS>

###### Abstract

Safely interacting with humans is a significant challenge for autonomous
driving. The performance of this interaction depends on machine learning-based
modules of an autopilot, such as perception, behavior prediction, and
planning. These modules require training datasets with high-quality labels and
a diverse range of realistic dynamic behaviors. Consequently, training such
modules to handle rare scenarios is difficult because they are, by definition,
rarely represented in real-world datasets. Hence, there is a practical need to
augment datasets with synthetic data covering these rare scenarios. In this
paper, we present a platform to _model_ dynamic and interactive scenarios,
_generate_ the scenarios in simulation with different modalities of labeled
sensor data, and _collect_ this information for data augmentation. To our
knowledge, this is the first integrated platform for these tasks specialized
to the autonomous driving domain.

##  1 Introduction

Dynamic interactions with humans is a significant challenge for autonomous
driving. Perception, behavior prediction, and planning modules of autonomous
vehicle software stack play crucial roles in interacting with humans on the
road. To train these machine learning-based modules, high quality labeled
sensor data with diverse and realistic dynamic behaviors is essential.
Handling rare scenarios for these modules is especially difficult because of
the, by definition, small proportion of such scenarios in real-world training
datasets. Hence, there is a practical need to augment these datasets with
synthetic data. However, in simulation, modeling dynamic and interactive
environments with complex spatiotemporal relations among multiple agents is
expensive, and is a serious barrier to scalable data augmentation.

We present a platform to model and generate dynamic scenarios, and collect
labeled sensor data with different sensor modalities. To our knowledge, this
is the first integrated platform for these tasks specialized to the autonomous
driving domain. We build on the Scenic scenario description language for
modeling and generating scenarios [1], extended to dynamic scenarios in [2].
Scenic does not provide flexible features for generating sensor data, and in
the context of data augmentation has only been used to generate static
2-dimensional RGB images [1, 3]. In this work, we add features to Scenic
enabling it to serve as a dynamic scenario data generation platform, with the
ability to: (1) configure arbitrary number and types of sensors on the ego
vehicle in simulation as they would be mounted on a real autonomous car, (2)
record data from these sensors along with corresponding ground-truth labels,
(3) extract specific types of sensor data from a recorded simulation, and (4)
visualize each type of recorded sensor data and labels. We open-source our
platform [4], along with 340GB of labeled sensor data generated from example
driving scenarios encoded in Scenic.

Figure 1: Overview of our platform’s pipeline.

##  2 Related Work

Several learning-based methods for scene generation have been proposed.
Hierarchical spatial priors between furniture are learned in [5, 6], which are
then integrated into hand-crafted cost functions for generating layouts of
indoor scenes. In [7], a generative model is proposed to sequentially add
objects into scenes. In [8], a dataset generator is parametrized with a neural
network. Generative Adversarial Networks (GANs) [9] have also used for data
augmentation [10, 11]. In contrast to such methods, Scenic does not require
training nor a pre-existing dataset to construct generative models. Also, as a
probabilistic programming language, Scenic can easily incorporate declarative
constraints, while producing data in an _explainable_ , programmatic fashion
requiring only a simulator.

Recently, domain-specific scenario description languages have been proposed
for autonomous cars. The Paracosm language [12] models dynamic scenarios with
a reactive and synchronous model of computation. Unlike Scenic, it lacks
probability distributions and declarative constraints, as well as constructs
like Scenic’s interrupts which enable easy customization of generic dynamic
behavior models. The Measurable Scenario Description Language (M-SDL) [13]
provides similar constructs as Scenic, but with some distinguishing features:
(i) Scenic provides a higher-level, declarative way to specify geometric
constraints, and (ii) Scenic is fundamentally a probabilistic programming
language, as opposed to M-SDL, whose distributions are optional.

##  3 Scenic: A Formal Dynamic Scenario Description Language

###  3.1 Scenario Modelling

Scenic [1, 2] is a formal scenario description language developed to model
static and dynamic, multi-agent scenarios to help designing and testing
autonomous systems such as self-driving cars. More precisely, Scenic is a
probabilistic programming language that enables users to specify distributions
over scenes (i.e., configurations of objects in space as well as their
features) and dynamic behaviors. Furthermore, users can construct objects in
an intuitive imperative style while simultaneously imposing declarative hard
(i.e., strict) and soft (i.e., probabilistic) constraints. Scenic’s concise,
readable syntax simplifies modelling _spatial_ and _temporal_ relationships.
This syntax provides building blocks to encode typical geometric relations
occurring in real-world driving scenarios, which would otherwise require
complex non-linear expressions and constraints, as well as temporal constructs
such as interrupts to model complex dynamic and reactive behaviors in a
modular fashion.

###  3.2 Scenario Generation

A Scenic program not only serves as a scenario description but also as a
source of synthetic data. An execution of the Scenic program samples an
environment according to the probabilistic distributions defined in the
program. At the beginning of program execution, Scenic samples the static
aspects of the environment, e.g. the initial positions and orientations of
objects, to define the initial condition for a simulation. As the simulation
runs, Scenic may sample further random variables pertaining to the dynamic
aspects of the environment, for example delays between events. The Scenic
program and an interfaced simulator communicate every simulation timestep to
generate reactive behaviors. The simulator sends all objects’ ground-truth
states (e.g., speed, heading, and position) to a server running the program.
Given this information, the server chooses and sends to the client (i.e., the
simulator) the next action (e.g. set throttle, steering angle, brake, etc) to
be simulated for each agent during the next simulation timestep.

For example, consider a scenario where a badly-parked car by the curb suddenly
pulls into the lane that the ego vehicle is driving on. A snippet of a Scenic
program modeling this scenario is shown in Figure 2. Data generated with this
program is shown in Figure 3.

Figure 2: Part of a Scenic program modeling a badly-parked car pulling into
the ego’s driving lane.

##  4 Platform Pipeline

The main contribution of this paper is a platform to generate any size of
synthetic dataset covering dynamic scenarios of interest. We used the CARLA
simulator [14] as a simulation and data generation engine for this platform.
However, Scenic is simulator-independent and can be interfaced with different
simulators to collect data if the they support relevant functionalities.

As illustrated in Figure 1, the input to the platform is a batch of Scenic
programs, each modeling a particular dynamic scenario. In addition, one can
configure any number of camera and lidar sensors, each with its own set of
parameters (e.g., sensor position, resolution and field of view for video
cameras, or rotational frequency for lidar). The duration of each simulation
as well as the number of simulations to run per scenario are adjustable. The
Scenic programs and configuration settings are compiled and sent to CARLA to
output a dataset of recorded scenarios, complete with ground-truth annotations
for each sensor. These are RGB video and lidar data with ground-truth
annotations for semantic segmentation, depth information, labeled point
clouds, and 2D/3D bounding boxes as shown in Figure 3. The details of how to
use our platform are available at [4].

Our implementation also comes with an API to browse the synthetic datasets
output by the platform, offering a way of programmatically interacting with
annotated data in order to train and evaluate machine learning models. For
instance, the API includes helper functions to project ground-truth 3D
bounding box annotations to 2D bounding boxes on any given camera (implicitly
computing the camera’s calibration matrix), and other functions to visualize
semantic point clouds with color labels.

##  5 Experiments

(a) Semantic Segmentation

(b) Lidar with Semantic Labels

(c) 2D Bounding Boxes

(d) Depth Image

(e) 3D Bounding Boxes

Figure 3: Collected Sensor Data and Labels using the Scenic program described
in Figure 2

To demonstrate the effectiveness of the platform, we created a synthetic
dataset of hundreds of simulated video recordings. This data were generated
from a variety of Scenic programs representing a range of different driving
scenarios. Many of these scenarios were adapted from the CARLA Autonomous
Driving Challenge Scenarios [15], which are derived from a pre-crash typology
reported by the National Highway Traffic Safety Administration [16]. Others
were chosen from the Voyage’s Open Autonomous Safety [17], an open-source
library containing commonly-encountered road scenarios used for AV safety
testing. Additional scenarios capturing high-density traffic situations, such
as bumper-to-bumper lines of cars and slow-moving intersectional left turns,
were also included. Each scenario was implemented as a Scenic program and
small variations such as background vehicles were added to the scenarios
provided in [15] and [17] in order to make the simulations more realistic.

For the synthetic dataset, an RGB camera was used to record front-facing
videos for the ego vehicle, along with ground-truth depth and semantic
segmentation data. In addition, a lidar sensor recorded 360-degree sweeps of
the scene, including semantic labels for each point. The simulator directly
reported 3D bounding boxes for objects in the scene, whose coordinates were
transformed to be relative to the ego. Figure 3 illustrates the provided
annotations on a single frame. The video camera and lidar sensor were placed
at the same location on the ego vehicle’s hood, exactly 2.4 meters off the
ground. The video camera recorded 15 frames per second at a resolution of
1280x720, with a 90 degree horizontal field of view (FoV). Per-pixel-depth
information and semantic labels were also captured. The lidar sensor ran at a
frequency of approximately 15 sweeps per second, recording about 22,000 points
per sweep. A complete description of all the semantic tags and exact sensor
configurations can be found at [4].

##  6 Conclusion

We presented a platform to model dynamic and interactive scenarios of
interest, generate such scenarios in simulation, collect annotated sensor
data, and visualize the recorded data. Our platform allows researchers and
developers to configure sensors as installed on their real autonomous car to
help them reduce the domain gap between the real world and the synthetic
dataset. We hope that our platform can assist in designing systems to safely
interact with humans on the road by allowing fine-grained control over the
distributions of different types of scenarios in the dataset.

## References

  * [1] D. Fremont, X. Yue, T. Dreossi, S. Ghosh, A. L. Sangiovanni-Vincentelli, and S. A. Seshia, “Scenic: A language for scenario specification and scene generation,” Programming Language Design and Implementation (PLDI), 2018\.
  * [2] D. J. Fremont, E. Kim, T. Dreossi, S. Ghosh, X. Yue, A. L. Sangiovanni-Vincentelli, and S. A. Seshia, “Scenic: A language for scenario specification and data generation.” http://arxiv.org/abs/2010.06580, 2020\.
  * [3] D. J. Fremont, J. Chiu, D. D. Margineantu, D. Osipychev, and S. A. Seshia, “Formal analysis and redesign of a neural network-based aircraft taxiing system with VerifAI,” in 32nd International Conference on Computer Aided Verification (CAV), July 2020.
  * [4] “Platform code base.” https://github.com/jayshenoy/scenic.
  * [5] K. Wang, Y.-A. Lin, B. Weissmann, M. Savva, A. X. Chang, and D. Ritchie, “Planit: Planning and instantiating indoor scenes with relation graph and spatial prior networks,” ACM Transactions on Graphics (TOG), vol. 38, no. 4, pp. 1–15, 2019.
  * [6] L. F. Yu, S. K. Yeung, C. K. Tang, D. Terzopoulos, T. F. Chan, and S. J. Osher, “Make it home: automatic optimization of furniture arrangement,” ACM Transactions on Graphics (TOG)-Proceedings of ACM SIGGRAPH 2011, v. 30,(4), July 2011, article no. 86, vol. 30, no. 4, 2011.
  * [7] S. A. Eslami, N. Heess, T. Weber, Y. Tassa, D. Szepesvari, G. E. Hinton, et al., “Attend, infer, repeat: Fast scene understanding with generative models,” in Advances in Neural Information Processing Systems, pp. 3225–3233, 2016.
  * [8] A. Kar, A. Prakash, M.-Y. Liu, E. Cameracci, J. Yuan, M. Rusiniak, D. Acuna, A. Torralba, and S. Fidler, “Meta-sim: Learning to generate synthetic datasets,” in Proceedings of the IEEE International Conference on Computer Vision, pp. 4551–4560, 2019.
  * [9] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” in Advances in neural information processing systems, pp. 2672–2680, 2014.
  * [10] X. Liang, Z. Hu, H. Zhang, C. Gan, and E. P. Xing, “Recurrent topic-transition gan for visual paragraph generation,” arXiv preprint arXiv:1703.07022, 2017\.
  * [11] M. Marchesi, “Megapixel size image creation using generative adversarial networks,” arXiv preprint arXiv:1706.00082, 2017.
  * [12] R. Majumdar, A. Mathur, M. Pirron, L. Stegner, and D. Zufferey, “Paracosm: A language and tool for testing autonomous driving systems,” 2019.
  * [13] Foretellix, “Measurable scenario description language.” https://www.foretellix.com/wp-content/uploads/2020/07/M-SDL_LRM_OS.pdf, 2020\.
  * [14] A. Dosovitskiy, G. Ros, F. Codevilla, A. M. López, and V. Koltun, “CARLA: an open urban driving simulator,” CoRR, vol. abs/1711.03938, 2017.
  * [15] “Carla autonomous driving challenge.” https://carlachallenge.org/.
  * [16] W. G. Najm, J. D. Smith, and M. Yanagisawa, “Pre-crash scenario typology for crash avoidance research,” 2007.
  * [17] “Voyage open autonomous safety.” https://oas.voyage.auto/.


# EPA Particulate Matter Data - Analyses using Local Control Strategy

Robert L. Obenchain and S. Stanley Young

(September 2022)

###### Abstract

Analyses of large observational datasets tend to be complicated and prone to
fault depending upon the variable selection, data cleaning and analytic
methods employed. Here, we discuss analyses of 2016 US environmental
epidemiology data and outline some potential implications of our Non-
parametric and Unsupervised Learning perspective. Readers are invited to
download the CSV files we have contributed to dryad and apply the analytic
approaches they think appropriate. We hope to encourage development of a
broad-based “consensus view” of the potential effects of Secondary Organic
Aerosols (Volatile Organic Compounds that have predominantly Biogenic or
Anthropogenic origin) within PM2.5 particulate matter on Circulatory and/or
Respiratory mortality. The analyses described here ultimately focus on the
question: “Is life in a region with relatively high air-borne Biogenic
particulate matter also relatively dangerous in terms of Circulatory and/or
Respiratory mortality?”

Keywords: Nonparametric Unsupervised Learning, Local Control Strategy,
Clustering as Matching, Permutation Distributions.

##  1 Introduction

When a researcher starts analyzing a large cross-sectional dataset downloaded
from the internet, it’s too late to worry about deficiencies in the
experimental design characteristics of observational data! It is what it is.
You could discard parts of the data to make the remainder “look better”: e.g.
look more blocked and/or balanced. But using less than essentially “all” of
the available data opens you up to accusations of “cherry-picking” or having
“gardened” the data to get the answer you want, Glaeser (2006).

One sound approach is to use R-functions and CRAN-packages that implement a
comprehensive and adaptive strategy that helps you to robustly “design” your
analyses of potentially diverse and/or confounded measurements, as illustrated
here using data from Pye et al. (2021). The approach illustrated here uses the
LocalControlStrategy R-package of Obenchain (2015-2019).

The scatter plot in Figure 1 displays $x=$ Bvoc model predictions
($\mu{g}/m^{3}$) from the EPA verses $y=$ AACRmort estimates from the CDC for
$2,980$ US Counties. The Bvoc component of PM2.5 consists of volatile organic
compounds that are primarily Biogenic (rather than Anthropogenic) in origin
within “Secondary Organic Aerosols”, Pye et al. (2021). AACRmort values are
“Age Adjusted Circulatory and/or Respiratory mortality” counts of deaths per
100,000 residents.

Figure 1:  The above scatter for 2980 US Counties within the contiguous 48
States and Washington DC displays considerable variability in $y=$ AACRmort
for any given value of $x=$ Bvoc. The five (red, blue or green) regression
quantiles fitted using the quantreg $R-$package, Koenker (2005-2022), display
considerable non-linearity. Thus, given only the Bvoc level of an otherwise
unspecified US County, it would clearly be unrealistic to expect a meaningful
point-prediction of its AACRmort count to result.

In natural products research, it is standard practice to examine fractions of
raw natural material to figure out where an “active agent” is located. There
typically are many fewer compounds in the active fraction, so chemical
identification becomes more feasible. Starting in 1993, evidence pointing to
particulate matter smaller than $2.5$ microns, thereafter known as PM2.5,
became the proposed active fraction in the 1952 London Fog death event
…according to (unpublished) work by British environmentalists.

In reality, PM2.5 is a complex mixture of over four thousand chemical
compounds which vary over both time and space (urban vs. rural, etc.) Pye et
al. (2021) appear to have followed a natural products strategy by arguing that
“Secondary Organic Aerosols” (SOAs) could be “the” active agents in PM2.5.
Unfortunately, in the almost 30 years since the first PM2.5 publications
sponsored by the US EPA, neither measures of PM2.5 itself, nor any fraction of
it, has been shown to be genuinely predictive of mortality at the level of
individual US Counties or Cities.

Young, Kindzierski and Randall (2021) suggest that the current EPA paradigm
that “PM2.5 is, by itself, lethal” is wrong. While the efforts of Pye et al.
(2021) may be noble, potential challenges to their proposals abound.

##  2 Predicted rather than Observed Measures of SOAs

Although actual measurements of SOA components in air samples are apparently
not currently feasible, various “models” for predicting their presence from
measurements of other pollutants have apparently been used since Volkamer et
al. (2006). The analyses presented in Pye et al. (2021) appear to be based on
EPA CMAQ models (2019), but Pye et al. (2021) do not discuss the assumptions
or limitations of those models.

For example, a potential issue is that CMAQ methods may be loosely based upon
Generalized Additive Model (GAM) concepts. While Pye et al. (2021) do
reference Wood (2003,2004,2022), their introduction states “We use multiple
linear regression…”. Since GAMs are a “composite” of individual models, some
parts could be fit in unspecified “non-standard” ways. All that we know as
outside observers is that their EPA models do a much better job of
“predicting” CDC measures of Circulatory and Respiratory mortality (raw or
age-adjusted) than we would have expected. In short, we would much rather have
actual “observed” measurements of secondary organic aerosols from validated
scientific instruments to contribute to dryad and analyze here.

###  2.1 Variables in the Data from the CDC, EPA & Other Federal Agencies

While there are $123$ variables (columns) and $2,980$ observations (rows) in
the larger CSV file that we have archived at dryad, Young and Obenchain
(August 2022), the analyses presented here actually use only the
$11-$variables described in Table 1 on page 4. These $11$ variables are in the
smaller CSV file in the same dryad archive that contains just $25$ variables
and $2,973$ observations without any missing values. Six US Counties had NA
codes on the rather important PREMdeath variable (FIPS = 8033, 30055, 31091,
46069, 46073 and 49031.) A $7^{th}$ NA code occurred for FIPS = 32011 in the
IncomIEQ variable, which turned out to be a relatively unimportant variable.

The analyses presented here will ultimately focus on “potential causes” of
Circulatory and Respiratory Mortality within individual US Counties. Our focus
on the AACRmort variable from the CDC (rather than CRmort) stems primarily
from our previous experience [Obenchain, Young and Krstic (2019)] where the
proportion of county residents over 65 was a key predictor of local mortality.
While we did perform some preliminary analyses using CRmort as the primary
$y-$Outcome, we found them to be both more complex and less self-consistent
than the corresponding AACRmort analyses that we present below.

Our somewhat surprising and unexpected findings concerning distinct components
of pmTOT (i.e. PM2.5) focus primarily upon Bvoc, although the Avoc and pmSO4
variables from the EPA were also used in our final analyses.

TABLE 1 – Variable Information
---
Name | Description | Range | Source
FIPS | Federal Information Processing code ($4$ or $5$ digits) | 1001-56045 | FED
CRmort | Crude Circulatory and/or Respiratory Mortality (per 100K) | 64.8-1564 | CDC
AACRmort | Age Adjusted Circulatory and/or Respiratory Mortality (per 100K) | 95.5-802.6 | CDC
pmTOT | Particles in Air with diameter less than 2.5 micrometres ($\mu$g/m3) | 2.06-14.32 | EPA
Bvoc | (predicted) Biogenic Volatile Organic Compounds in pmTOT | 0.26-3.31 | EPA
Avoc | (predicted) Anthropogenic Volatile Organic Compounds in pmTOT | 0.23-2.89 | EPA
pmSO4 | Sulfate Compounds in pmTOT | 0.39-1.62 | EPA
ASmoke | Adult Smoking Fraction | 0.007-0.412 | EPA
PREMdeath | Premature Death Rate | 0.03-0.66 | EPA
ChildPOV | Children Living in Poverty (per 100K) | 2853-36469 | EPA
IncomIEQ | Income In-Equality Rating | 2.9-8.9 | EPA

Although we ultimately focus on Local (within Cluster) models that are highly-
flexible and statistically robust, an initial analysis of a simple linear
(regression) model was used to assess the “overall” extent of ill-conditioning
(confounding) among variables; see Figure 2 on page 5. This Trace display
shows fitted regression coefficients for the Efficient Generalized Ridge
Regression shrinkage path of Obenchain (2022). This linear model attempts to
predict $AACRmort$ using $Bvoc,Avoc,pmSO4,PREMdeath,ASmoke,ChildPOV$ and
$IncomIEQ$. While the order in which these seven “right-hand-side” variables
are specified is irrelevant in classical linear regression, the first three
are EPA environmental variables while the last four quantify socioeconomic
characteristics of county residents.

Figure 2:  Following a modest initial extent of shrinkage to $m=0.51$, MSE-
risk optimal relative-magnitudes for $\beta-$coefficient estimates are
determined. In other words, overall ill-conditioning (redundancy) among these
seven $X-$variables appears minor. Since our main LC Strategy analyses will
use AACRmort as the $y-$Outcome measure, we will focus on Bvoc as our primary
(potentially causal) Exposure variable. After all, the EPA Avoc and IncomIEQ
variables have only small and negative coefficient estimates. The genuinely
“bad news” from this classical linear model is that its overall
$R^{2}-$Coefficient of Determination is rather low (coincidentally, also
$0.51$). The “local” models developed here using LC Strategy will be non-
linear, distribution-free and will provide much better fits.

##  3 Applying Local Control Strategy

While LC Strategy can be used by single researchers working alone, it may also
be used by groups of collaborating researchers. In fact, LC Strategy may be an
ideal approach for use when the researches involved in a collaboration have
“odd combinations” of qualifications …such as: diverse experiences and unique
perspectives on the overall issues involved. Use of LC Strategy plus access to
data then makes it possible to seriously address the question: Can any shared,
consensus position be reached?

Our discussion of LC Strategy below will use the terminology introduced in
Table 2, also on page 5.

TABLE 2 – Basic LC Strategy Terminology
---
Term | Description
Experimental Unit | An individual US County or Parish.
Cluster | A subgroup of Units that are relatively well-matched on their given $X-$characteristics.
Cycle | One “pass” through the data involving at least the Aggregate and Confirm phases of analysis.
Aggregate Phase | The Phase of LC Strategy where new sets of Clusters are formed.
Confirm Phase | The Phase where a given Aggregation is shown to be either Ignorable or to provide
| Added Information.
Explore Phase | The Phase where Aggregations are ranked on Added Information and a decision of
| either Continue or Stop is made.
Reveal Phase | The final Phase of LC Strategy where the Variance-Bias Trade-Offs associated with
| the best Aggregation are Revealed and Discussed.

###  3.1 Forming the Hierarchical Clustering Tree

The first step in LC Strategy is to form a hierarchical clustering dendrogram.
The LocalControlStrategy::LCcluster() function calls stats::prcomp() to
calculate Mahalanobis distances using standardized and orthogonal Principal
Coordinates, Obenchain (1971), Rubin (1980). These coordinates optimally
quantify “dissimilarity” among US Counties using their possibly confounded
$X-$characteristics. A dendrogram (tree-like structure) can usually be “cut”
by a horizontal line to produce any desired number of individual clusters. Our
dendrogram, displayed in Figure 3, uses the default “ward.D” clustering
algorithm.

Figure 3:  This Hierarchical Clustering Dendrogram can be “cut” by a
horizontal line at any strictly positive Height. For example, note that the
displayed red line produces 10 clusters.

Either the divisive cluster::diana() method or any of six agglomerative
alternatives (ward.D2, complete, average, mcquitty, median, or centroid) could
have been specified. However, “single-linkage” clustering is not appropriate
for use in LC Strategy because it can produce clusters that are more like
“strings” than like compact subgroups of most similar experimental units.

###  3.2 Confirming that a given Number of Clusters provide “Added
Information”

Every new Aggregate Phase of LC Strategy is followed immediately by a Confirm
Phase analysis asking whether the new “LRC distribution” could possibly be
generated in some Completely Random way. Actually, it is easy to generate
purely random LRC distributions that use both the same number of clusters and
the same cluster sizes as any given “new” LRC distribution.

Figure 4 displays the numerical example with $K=50$ Clusters that we will use
to illustrate confirm phase testing. We are “jumping the gun” here in the
sense that we will argue next (Subsection $3.3$) that $K=50$ is also the
“optimal” number of Clusters to use in applying LC Strategy.

LRC distributions are always discreet, while the traditional test that uses
the Kolmogorov-Smirnov D-statistic is only valid for comparing two absolutely
continuous distributions. Thus, LC Strategy uses a (computationally intensive)
permutation test, Welch (1990), that is non-parametric. Specifically, the LC
Strategy confirm() function easily simulates the gray eCDF displayed in Figure
4 and computes the resulting Kolmogorov-Smirnov D-statistic. These
calculations take less than $5$ seconds.

Note that the gray NULL eCDF appears to be rather smooth in Figure 4. After
all, it results from pooling $100$ sets of LRC estimates from $50$ pseudo-
clusters that were formed in a purely random and meaningless way. Thus the
gray eCDF contains at most $5,000$ quite small steps. To truly provide “added
information”, the distribution underlying the Blue eCDF from clusters of “well
matched” US Counties must be clearly different from the purely random
distribution (from random pseudo-clusters) represented by the gray eCDF.

The primary objective of the KSperm() function in LC Strategy is to simulate
an appropriate p-value, “adjusted” for ties, when testing the NULL hypothesis
that all x-Covariates used to form clusters are essentially Ignorable. These
computations, using $1,000$ replications, take roughly $77$ seconds on a
desktop PC.

Figure 5 displays the NULL distribution where the largest observed NULL
D-statistic is only $0.237$. This simulated NULL D-statistic is dwarfed by the
Observed $D-$statistic of $0.7903$ in Figure 4. In other words, the Blue
empirical CDF variable, named “LRC50”, provides abundant Added Information in
the sense that it is clearly “shifted to the left” in Figure 4 relative to the
simulated gray NULL eCDF. The Blue eCDF is thus Not Ignorable! Furthermore,
note that some Local Rank Correlation estimates from $K=50$ clusters are
clearly negative.

Figure 4:  This plot displays a pair of empirical Cumulative Distribution
Functions (eCDFs) and addresses the question of whether they are statistically
different. Note that the Kolmogorov-Smirnov D-statistic of $0.7903$ (maximum
vertical separation between eCDFs) occurs here at $0.30$ on the horizontal
axis …marked by the dashed red line. Figure 5:  Since the traditional
Kolmogorov-Smirnov test compares eCDFs from a pair of Continuous
Distributions, it is severely biased when used to compare Discrete
Distributions. Thus, simulation is used in LC Strategy to make valid
comparisons between eCDF estimates. Since the largest observed NULL
D-statistic was only $0.237$, the Observed D of $0.7903$ almost surely has an
actual $p-$Value of much less than $0.001$ (one in a thousand).

###  3.3 Choosing the Number of Clusters to Use in LC Strategy

Especially when cross-sectional data are observational, the data typically
contain isolated “local” effects as well as overall (global) confounding from
potentially highly correlated variables. The Nonparametric and Unsupervised
“pre-processing” methods that characterize LC Strategy address these issues.
Traditional methods based on multiple regression are typically Global,
Parametric and Supervised. They tend to make much stronger (and possibly
unrealistic) assumptions that yield questionable global predictions and
potential extrapolations.

In sharp contrast, LC Strategy summarizes findings from several separate and
flexible “local” models that are fit within clusters of “most similar” and/or
relatively “well-matched” experimental units. The number of such local models
and separate clusters of experimental units being used in a given Cycle is
denoted by the symbol K.

LC Strategy embraces many of the basic concepts and experiences outlined in
Stuart (2010), van der Laan and Rose (2010) and Stang et al. (2010). Published
analyses using LC Strategy include: Obenchain and Young (2013); Lopiano,
Obenchain and Young(2014); Young, Smith and Lopiano (2017) and Obenchain,
Young and Krstic (2019).

The LC approach consists of Four Phases of analysis called, respectively :
Aggregate, Confirm, Explore and Reveal. Each analysis “cycle” through the
available data typically uses a different Aggregation of experimental units
into a larger and larger number, $K$, of smaller and smaller “Clusters”.
Potential results from each new cycle are then either discarded or else
Confirmed to provide Added Information …in the sense of being clearly
different from “purely random” in ways discussed below.

The implicit statistical model used in LC Strategy is simply one-way ANOVA.
The Within-Block measures of effect-size when the $e-$variable is an Exposure
(with many more than just two levels) are “Local Rank Correlations” (LRCs)
between Exposures and $y-$Outcomes. A single cluster containing $N$ US
Counties thus contributes $N$ identical LRC estimates to the overall (across
cluster) distribution of LRC “effect-sizes”.

Figure 6:  Here we see $7$ I-plots of the overall distribution of LRC effect-
sizes as the number of clusters, $K$, being formed increases from $1$ to
$200$. Note that there is clear evidence that Variance$-$Bias Trade$-$Offs are
occurring for $K\geq 50$. For example, the Vertical inter-quartile ranges of
local distributions then steadily increase, and the median of local
effect$-$sizes is minimized at $K=100$. However, the vertical spread between
upper and lower extremes of LRC distributions start becoming excessive for
$K>50$. Our recommended “compromise” choice is thus $K=50$.

Ultimately, we chose to use $K=50$ Clusters of US Counties formed using only
$6-$Key $X-$confounder variables, out of the more than $100$ that were
available. These $50$ clusters represent mutually exclusive and exhaustive
statistical “Blocks” of relatively well-matched experimental units, and our
resulting “best” Effect-Size distribution quantifies a potentially causal
relationship between the numeric $y-$Outcome variable AACRmort and the
$e-$Exposure variable Bvoc.

Our choice of $K=50$ clusters was determined within the Explore phase when
$Box-Plots$ showed that the distribution of estimated “effect sizes” would
become excessively variable for $K>50$. Specifically, any analyst applying LC
Strategy would examine the “LCcompare()” graphic displayed in Figure 6 on page
10 and choose $K=50$ for the very reasons outlined in the caption of that
figure.

Figure 7:  This histogram shows the number of US Counties within each of the
50 Clusters formed using observational data from the EPA, CDC and US Census.
In a prospective study (e.g. a well “designed” experiment), all “blocks” of
well-matched experimental units are usually of nearly equal size. But that’s
an unlikely “ideal” situation when analyses must rely on observational data.
Figure 8:  Here we see that the distribution of Local Rank Correlations (LRCs)
within our 50 Clusters spans the range from $-0.30$ to $+0.60$. Since this
distribution contains only $50$ distinct numerical values, it certainly does
not (and cannot) look much like a “continuous” parametric distribution. Figure
9:  Here are plots of Bvoc versus AACRmort within four of the $50$ Clusters.
The straight blue lines are least-squares fits, while the red curves are
“loess” smooths (span = $0.75$) for predicting AACRmort from Bvoc. While the
(Spearman) LRCs within the two left-hand Clusters ($11$ and $40$) are
positive, both right-hand clusters ($26$ and $49$) are typical of the negative
LRC estimates observed for a total of $421$ US Counties ($21.3\%$ of the
$2,973$ US Counties without missing data).

##  4 LC “Reveal” Stage: randomForests, Partial Dependence Plots and a Single
RP Tree

We can now stop our “Explore” phase analyses and enter a (final) “Reveal”
phase. Specifically, we now look across, rather than only within, the $K=50$
clusters of relatively well-matched US Counties that provide the most
(potentially causal) Added Information about relationships between the
AACRmort and Bvoc “measures” of mortality and air pollution, respectively.
This new information is embedded within the LRC50 variable that takes on only
one of $50$ distinct numerical values between $-0.699$ and $+0.600$ for each
US County.

###  4.1 Interaction Effects Abound

Figure 9 contains plots of Bvoc versus AACRmort for $4$ of our $50$ Clusters.
The two left-hand clusters have positive LRC-estimates and blue least-squares
lines that slope upwards, while the two right-hand clusters have negative LRC-
estimates and blue least-squares lines that slope downwards. Actually, $421$
US Counties ($21.3\%$) of all $2,973$ US Counties without missing data fall
into the $10$ out of $50$ Clusters (i.e. $20\%$) that yield negative LRC-
estimates. In other words, almost $80\%$ of US Counties are in the $80\%$ of
Clusters with strictly positive LRC-estimates. After all, the Rank-Correlation
for the overall Bvoc vs. AACRmort scatter displayed in Figure 1 is $+0.474$.
This made us wonder whether larger clusters might tend to produce positive LRC
estimates?

An extreme illustration of this is displayed in the lower-right panel of
Figure 10. Data from the $237$ Counties within the $6$ clusters with the most
Negative LRCs are merged together there, yielding a scatter with clearly
positive slope (and an LRC of $+0.267$.) Meanwhile, each of the other $3$
panels shows scatters for $2$ out of the $6$ clusters with most negative LRCs.
Clearly, clustering of US Counties on their confounded $X-$characteristics
predictive of AACRmort truly does matter!

Figure 10:  When AACRmort and Bvoc data from the 6 clusters with most Negative
LRC estimates are merged together, much of their cluster “identity” gets lost.
This is clear in the lower-right panel: their “merged trend” is then clearly
Positive (shown in Yellow). Each of the other three panels show AACRmort vs.
Bvoc scatters for two clusters with negative LRCs using a total of six
different colors.

###  4.2 Insights from a randomForest of $500$ Trees

The highly detailed captions for the next eight Figures (11 through 18)
outline our interpretations of the relationships revealed by Partial
Dependence Plots (PDPs) within our “randomForest” of $500$ tree models.

Figure 11:  The single most important predictor of LRCs is ASmoke
($\%$IncMSE$=71.28$). Adult Smoking rates within the $.15$ to $.20$ range are
most predictive of positive LRCs between AACRmort and Bvoc. For some unknown
reason, smoking rates above $.25$ correspond to somewhat lower but still
relatively high LRCs. Figure 12:  The second most important predictor of LRCs
is ChildPOV ($\%$IncMSE$=59.46$). When Children in Poverty rates are less than
$.30$, LRCs are strongly positive. But LRCs then drop and can be negative when
ChildPOV $\geq.50$. In fact, ChildPOV is the only predictor associated with
strictly negative LRC predictions …where AACRmort counts tend to decrease as
Bvoc levels increase. Figure 13:  The third most important predictor of LRCs
is PREMdeath ($\%$IncMSE$=56.31$). PREMdeath is the (raw) count of premature
deaths per $100,000$ county residents. As this rate increases from $5,000$ to
$15,000$, LRCs steadily increase from $0.15$ to about $0.28$. LRCs then
approach $0.30$ for all PREMdeath counts exceeding $20,000$. A high PREMdeath
count appears to be a truly insidious indicator of poor Circulatory and/or
Respiratory Health. Figure 14:  The fourth most important predictor of LRCs
(Within-Cluster Spearman rank correlations between AACRmort and Bvoc) is Bvoc
itself ($\%$IncMSE$=46.40$). As local Bvoc levels increase from $0.2$ to
$1.5$, LRCs tend to steadily increase from $0.15$ to a maximum of $0.24$ and
remain at this highest level for further increases in Bvoc level. Figure 15:
The fifth important predictor of LRCs is pmSO4 ($\%$IncMSE$=37.82$). This
PDPlot suggests, curiously, that, while LRCs steadily increase over the range
$0.4-1.4$, they sharply drop when pmSO4 exceeds $1.4$. Figure 16:  The sixth
important predictor of LRCs is Avoc ($\%$IncMSE$=25.26$). Here, LRCs tend to
mostly increase from about $0.12$ to $0.20$ over the initial range of $0.1\leq
Avoc\leq 1.5$, then remain above $0.22$ when Avoc exceeds $1.8$. The general
shape of this PDP is somewhat like that for the Bvoc predictor, Figure 14.
Figure 17:  A seventh predictor of LRCs (between AACRmort and Bvoc) is
AACRmort itself: $\%$IncMSE$=17.79$. Since the vertical range in this plot is
quite small ($0.014$), any suggestion of “curvature” here is relatively
unimportant. Figure 18:  The last (and least important) predictor of within-
cluster LRCs is IncomIEQ $\%$IncMSE$=16.15$. LRCs tend to be near to $0.185$
for $3\leq IncomIEQ\leq 5.5$, then steadily drop to $0.165$ for IncomIEQ above
$8$.

Table 3 summarizes the key characteristics of our eight PDP plots. The top 4
$X-$confounder variables predictive of AACRmort include the three variables
inducing highest Incremental Node Purity.

TABLE 3 – Importance Statistics for $8$ Variables
---
Variable | %IncMSE | IncNodePurity
ASmoke | 71.28070 | 16.611216
ChildPOV | 59.45606 | 17.160771
PREMdeath | 56.31122 | 11.765453
Bvoc | 46.40073 | 16.198999
pmSO4 | 37.81722 | 12.071959
Avoc | 25.25779 | 14.658269
AACRmort | 17.79182 | 4.878167
IncomIEQ | 16.14760 | 5.173446

###  4.3 Insights from Model-Based Recursive Partitioning

We used the “model based” algorithms of Zeileis, Hothorn and Hornik (2008) to
develop the single Recursive Partitioning tree displayed in Figure 19. Note
that this approach searches for “optimal intervals” of Bvoc values within
which observed LRCs can be predicted from local ASmoke levels using simple
linear regression. Summary statistics for the $7$ Bvoc intervals that resulted
are listed in Table $4$ on page 18.

Figure 19:  This figure suggests that high ASmoke rates are particularly
harmful within the two lowest Bvoc intervals. Within four of the five higher
Bvoc intervals towards the right, LRCs (between AACRmort and Bvoc) are
surprisingly predicted to decrease as local ASmoke rates increase. TABLE 4 –
Statistics for LRC Prediction within $7$ Bvoc Intervals
---
Final | LRC | ASmoke | Bvoc | Bvoc | US
Node | Intercept | Slope | Range | Level | Counties
4 | -0.00411 | +1.09902 |  $\leq$ 0.313 | Low | 126
5 | -0.13183 | +1.42106 | 0.313-0.439 | Low | 379
=== | === | === | === | === | ===
7 | 0.33623 | -1.24350 | 0.439-0.805 | Medium | 862
8 | 0.12447 | +0.29020 | 0.805-1.010 | Medium | 323
=== | === | === | === | === | ===
10 | 0.45914 | -1.01360 | 1.010-1.718 | High | 750
12 | 0.86823 | -3.21152 | 1.718-2.237 | High | 311
13 | 1.31711 | -5.09834 |  $>$ 2.237 | High | 222

##  5 Summary

Personal computers have helped shape statistical theory as well as its
practice over the last 50 years. Freely available software can provide
computational and visual fast-tracks into the strengths and weaknesses of
alternative statistical methods. For example, the lower right-hand panel of
Figure 10 showed just how misleading a scatter-plot can be that ignores
cluster-membership. Our clusters contain US Counties with most similar
$X-$characteristics, so they compare only “apples” with “apples”, “potatoes”
with “potatoes”, etc, etc. This strategy leads to either fair comparisons when
the exposure variable is binary [Lopiano, Obenchain and Young (2014)], or to
cogent local inferences when that exposure (Bvoc) is both continuous and
correlated with potential confounding variables (Avoc, pmSO4, ChildPOV, etc,
etc.)

The eight PDP plots that are displayed in Figures 10 to 17 and summarized in
Table 3 on page 14 have encouraged us to focus our attention primarily on the
top-four predictors of LRCs; Bvoc is number four. Unfortunately, neither
$pmSOA\equiv Avoc+Bvoc$ nor either of its two components can apparently be
actually “measured” today …using existing scientific instruments with
reasonable accuracy.

Bvocs include Terpenes from trees and grass, and EPA Bvoc predictions are, at
least on average, positively correlated with corresponding local CDC CRmort
counts, before and after being Age Adjusted. How will US environmentalists
react to this “news”? Luckily, our single RP tree in Figure 19 does suggest
that levels of Bvoc $>1.01$ can be protective against high AACRmort when
ASmoke rates are above roughly $20\%$.

The future of NU-Learning approaches for analysis of cross-sectional
observational data strikes us as being quite promising. Unfortunately,
traditional (Parametric and Supervised Learning) models used in Propensity
Score estimation tend to be global (rather than local) and may make strong and
potentially unrealistic assumptions.

Most interesting questions don’t have an answer that is clearly “Black” or
“White”. Data from well-designed and well-executed experiments can illuminate
some key issues. Meanwhile, large collections of observational data can reveal
interesting real-world “shades of gray”.

Figure 20:  This US map, Kahle and Wickham (2013), shows variation in EPA Bvoc
predictions across US Counties using six shades of Green. Low Bvoc levels
receiving the lightest shade of Green dominate the US “Great Plains” and Rocky
Mountain regions where agricultural production is either low or supported by
irrigation. Note that a vast majority of counties receive the same Bvoc color
as their adjacent counties. Note also that 151 US Counties fall into the
$7^{th}$ “NA” category and receive the see-through “background” color. Figure
21:  This map shows variation in AACRmort counts per $100,000$ residents using
six shades of Red. These adjusted counts are certainly not smooth functions of
geographic latitude and longitude; variation in AACRmort counts may result
primarily from variation in the socioeconomic characteristics of local
residents. After all, the vast majority of adjacent counties receive a
somewhat different or very different shade of Red. Figure 22:  On this map
showing variation in Local Rank Correlations (LRCs) across US Counties, six
shades of Blue identify LRC levels of association between AACRmort and Bvoc
within our 50 clusters of US counties well-matched in the $7-$dimensional
space of potential confounded characteristics. Note that all adjacent counties
receiving different shades of Blue must be in different Clusters and thus have
somewhat different characteristics. All $421$ counties with negative LRC
estimates receive the lightest shade of Blue. Again, 151 counties without data
receive the see-through “background” color.

Our use of LC Strategy in our analysis of rather voluminous and detailed data
highlights a truly key advantage of our computer intensive approach. Namely,
this approach generates many simple yet highly informative graphical displays
that range from simplistic to quite detailed. Using a wide variety of displays
can guide further analyses in unanticipated directions. For example, note the
great detail embedded within Figures 19 through 22. Using a PDF reader with
Adobe zoom capability enables comparisons among US counties that are “similar”
in terms of either Cluster Membership in seven dimensions (Bvoc, Avoc, pmSO4,
PREMdeath, ASmoke, ChildPOV and IncomIEQ) or Geography in two dimensions (lat
and long).

Our main contributions and forward-looking recommendations are:

(1.) We have made CSV files publicly available that contain data as much
“like” the CDC and EPS data used by Pye et al. (2021) as is possible by
private US citizens.

(2.) Use of LC Strategy in our analyses of these data has yielded findings
that are fundamentally different and quite possibly more relevant and
technically sound than those of Pye et al. (2021). Since interaction effects
abound within the data, it’s not surprising that local models provide better
fits than “smoothly” continuous global models.

(3.) It became clear in our earliest data analyses that Bvoc would be a better
predictor of AACRmort than Avoc or $pmSOA\equiv Avoc+Bvoc$. For example, the
(Pearson) correlation $cor(AACRmort,Bvoc)=0.4589$ is larger than both
$cor(AACRmort,Avoc)=0.2489$ and $cor(AACRmort,pmSOA)=0.4253$. Furthermore, all
Avoc $\beta-$coefficient estimates are negative in Figure 2. Did Pye et al.
(2021) avoid talking about Bvoc effects primarily because EPA regulation of
Bvoc would have no clear organization to hold responsible?

(4.) What, exactly, is the EPA explanation of why their CMAQ models (2019)
appear to be much better predictors of AACRmort (or CRmort before Age
Adjustment) than models from academic research traditionally funded by the
EPA? For example, does CMAQ use information on variation in socioeconomic
characteristics of residents? How and when will pmSOA and its components
become accurately measurable? We recommend considerable skepticism while
waiting for answers to these sorts of questions!

### R-code and Data Availability

A file named LCSonEPA$\\_$Code.pdf can be freely downloaded from
http://localcontrolstatistics.org. This PDF contains all R-code needed to
perform the analyses and display example figures. A variety of introductory
information about LC Strategy can also be freely downloaded from this internet
site.

Our reconstruction of all 123 variables used by Pye et al. (2021) can be
downloaded from https://datadryad.org. We merged files from the CDC and EPA
using FIPS codes and created CSV files with shorter variable names. The larger
“CDC$\\_$EPA$\\_$2022.csv” file contains $123$ variables while
“AnalysisFile.csv” contains only the $25$ variables either considered for or
actually used in the analyses presented in this paper. The smaller file
contains data on $2,973$ US Counties without “missing values”, while the
larger file contains everything in the smaller file, except the Cluster
membership number and LRC value, for $2,980$ US Counties.

### Conflict of Interest

The authors state no conflicts of interest.

## References

  * [1] Breiman, L. (2001), “Random Forests”, Machine Learning, 45, 5$-$32\. https://doi.org/10.1023/A:1010933404324.
  * [2] Breiman, L. (2002), “Manual On Setting Up, Using, And Understanding Random Forests, V3.1”. https://www.stat.berkeley.edu/~breiman/Using_random_forests_V3.1.pdf.
  * [3] Centers for Disease Control and Prevention National Center for Health Statistics. Compressed Mortality File 1999-2016 on CDC WONDER Online Database, released June 2017. Data are from the Compressed Mortality File 1999-2016 Series 20 No.2U, 2016, as compiled from data provided by the 7 vital statistics jurisdictions through the Vital Statistics Cooperative Program. Accessed by Obenchain on 19 April 2022, http://wonder.cdc.gov/cmf-icd10.html.
  * [4] Friedman, J. (2001), “Greedy function approximation: the gradient boosting machine”, Annals of Statistics, 29(5), 1189$-$1232, https://doi.org/10.1214/aos/1013203451.
  * [5] Glaeser, E.L. (2006), “Researcher incentives and empirical methods”, https://www.nber.org/system/files/working_papers/t0329/t0329.pdf.
  * [6] Hothorn, T. and Zeileis, A. (2015), “partykit: A Modular Toolkit for Recursive Partytioning in R.”, Journal of Machine Learning Research, 16, 3905$-$3909\. https://jmlr.org/papers/v16/hothorn15a.html.
  * [7] Hothorn, T., Seibold, H. and Zeileis, A. (2015-2022), “partykit: A Modular Toolkit for Recursive Partytioning”, ver 1.2-16, https://CRAN.R-project.org/package=partykit.
  * [8] Kahle, D. and Wickham, H. (2013), “ggmap: Spatial Visualization with ggplot2.”, The R Journal, 5, 144$-$161\. http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf.
  * [9] Koenker, R. (2005 - 2022), “quantreg: Quantile Regression”, ver 5.94, https://CRAN.R-project.org/package=quantreg.
  * [10] Liaw, A. and Wiener, M. (2002 - 2022) “randomForest: Breiman and Cutler’s Random Forests for Classification and Regression”, ver 4.7-1, https://CRAN.R-project.org/package=randomForest.
  * [11] Lopiano, K.K., Obenchain, R.L. and Young, S.S. (2014), “Fair Treatment Comparisons in Observational Research”, Statistical Analysis and Data Mining, 7, 376$-$384\. https://doi.org/10.1002/sam.11235.
  * [12] Obenchain, R.L. (1971), “Multivariate Procedures Invariant under Linear Transformations, Annals of Mathematical Statistics, 42: 1569$-$1578\. https://doi.org/10.1214/aoms/1177693155.
  * [13] Obenchain, R.L. (2005-2022), “RXshrink: Maximum Likelihood Shrinkage using Generalized Ridge or Least Angle Regression Methods”, ver 2.1, https://CRAN.R-project.org/package=RXshrink.
  * [14] Obenchain, R.L. (2015-2019), “LocalControlStrategy: Local Control Strategy for Robust Analysis of Cross-Sectional Data”, ver 1.3.3, https://CRAN.R-project.org/package=LocalControlStrategy.
  * [15] Obenchain, R.L. (2022), “Efficient Generalized Ridge Regression, Open Statistics, 3: 1$-$18\. https://doi.org/10.1515/stat-2022-0108.
  * [16] Obenchain, R.L. and Young, S.S. (2013), “Advancing statistical thinking in health care research”, Journal of Statistical Theory and Practice, 7, 456$-$469\. https://doi.org/10.1080/15598608.2013.772821.
  * [17] Obenchain, R.L., Young, S.S. and Krstic, G. (2019), “Low-level Radon Exposure and Lung Cancer Mortality”, Regulatory Toxicology and Pharmacology, 107, https://doi.org/10.1016/j.yrtph.2019.104418.
  * [18] Pye, H.O.T., Ward-Caviness, C.K., Murphy, B.N., Appel, K.W. and Seltzer, K.M. (2021), “Secondary organic aerosol association with cardiorespiratory disease mortality in the United States”. Nature Communications 12, 7215. https://doi.org/10.1038/s41467-021-27484-1.
  * [19] R Core Team. (2022), “R: A language and environment for statistical computing”. https://www.R-project.org/
  * [20] Rubin, D.B. (1980), “Bias reduction using Mahalanobis metric matching”. Biometrics, 36, 293$-$298\. https://www.jstor.org/stable/2529981.
  * [21] Stang, P.E., Ryan, P.B., Racoosin, J.A., Overhage, J.M., Hartzema, A.G., Reich, C., Welebob, E., Scarnecchia, T. and Woodcock, J. (2010), “Advancing the Science for Active Surveillance: Rationale and Design for the Observational Medical Outcomes Partnership”. Annals of Internal Medicine, 153, 600$-$606\. https://doi.org/10.7326/0003-4819-153-9-201011020-00010.
  * [22] Stuart, E.A. (2010), “Matching methods for causal inference: A review and a look forward”. Statistical Science, 25, 1$-$21\. https://doi:10.1214/09-STS313.
  * [23] U.S. EPA Office of Research and Development. (2019), “CMAQ (Version 5.3.1)”, https://doi.org/10.5281/zenodo.3585898.
  * [24] Volkamer, R., Jimenez, J.L., San Martini, F., Dzepina, K., Qi, Z., Salcedo, D., Molina, L.T., Worsnop, D.R. and Molina, M.J. (2006), “Secondary organic aerosol formation from anthropogenic air pollution: Rapid and higher than expected”. Geophys. Res. Lett. https://doi.org/10.1029/2006GL026899.
  * [25] Welch, W.J. (1990), “Construction of permutation tests.” Journal American Statistical Association, 85, 693$-$698\. https://doi.org/10.1080/01621459.1990.10474929.
  * [26] Wickham, H. (2016), “ggplot2: Elegant Graphics for Data Analysis”. New York: Springer-Verlag.
  * [27] Wood, S. N. (2003), “Thin plate regression splines”. J. R. Stat. Soc. B. Met. 65, 95$-$114\. https://doi.org/10.1111/1467-9868.00374.
  * [28] Wood, S. N. (2004), “Stable and efficient multiple smoothing parameter estimation for generalized additive models”. Journal American Statistical Association, 99, 673$-$686\. https://doi.org/10.1198/016214504000000980
  * [29] Wood, S. N. (2022), “mgcv: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation”, https://CRAN.R-project.org/package=mgcv.
  * [30] van der Laan, M. and Rose, S. (2010), “Statistics ready for a revolution: Next generation of statisticians must build tools for massive data sets”. AMStat News, September, 38–39. https://magazine.amstat.org/blog/2010/09/01/statrevolution/.
  * [31] Young, S.S., Kindzierski, W. and Randall, D. (May 2021), “Shifting Sands: Unsound Science and Unsafe Regulation. Report $\\#1$: Keeping Count of Government Science – P-Value Plotting, P-Hacking, and PM2.5 Regulation”. National Association of Scholars. https://files.eric.ed.gov/fulltext/ED616199.pdf.
  * [32] Young, S.S., Obenchain, R.L. (August 2022), “EPA Particulate Matter Data”. https://doi.org/10.5061/dryad.63xsj3v58.
  * [33] Young, S.S., Smith, R.L. and Lopiano, K.K. (2017), “Air Quality and Acute Deaths in California, 2000$-$2012”, Regulatory Toxicology and Pharmacology, 88, 173$-$184\. https://doi.org/10.1016/j.yrtph.2017.06.003.
  * [34] Zeileis, A., Hothorn, T. and Hornik, K. (2008). “Model-Based Recursive Partitioning.” Journal of Computational and Graphical Statistics, 17, 492$-$514\. https://doi.org/10.1198/106186008X319331.


# Calibration of luminosity correlations of gamma-ray bursts using quasars

Sarveshkumar Purohit<EMAIL_ADDRESS>Shantanu Desai
<EMAIL_ADDRESS>Department of Artificial Intelligence, Indian Institute of
Technology, Hyderabad, Telangana-502284, India  Department of Physics, Indian
Institute of Technology, Hyderabad, Telangana-502284, India

###### Abstract

In order to test the efficacy of Gamma-ray Bursts (GRBs) as cosmological
probes, we the characterize the scatter in the correlations between six pairs
of GRB related observables, which have previously also been studied in Tang et
al. (2021). However, some of these observables depend on the luminosity
distance, for which one needs to assume an underlying cosmological model. In
order to circumvent this circularity problem, we use X-ray and UV fluxes of
quasars as distance anchors to calculate the luminosity distance in a model-
independent manner, which in turn gets used to calculate the GRB-related
quantities. We find that all the six pairs of regression relations show a high
intrinsic scatter for both the low and high redshift sample. This implies that
these GRB observables cannot be used as model-independent cosmological probes.

##  I Introduction

Gamma-ray bursts (GRBs) are single-shot explosions located at cosmological
distances, which were first detected in 1960s and have been observed over ten
decades in energies from keV to over 10 TeV range Kumar and Zhang (2015); Yu
et al. (2022). They are located at cosmological distances, although a distinct
time-dilation signature in the light curves is yet to be demonstrated Singh
and Desai (2022). Because of their high energies and cosmological distances,
they also have proved to be very good probes of fundamental physics, such as
testing of Lorentz invariance violation and quantum gravity Desai (2023). GRBs
are traditionally divided into two categories based on their durations, with
long (short) GRBs lasting more (less) than two seconds Kouveliotou et al.
(1993). Long GRBs are usually associated with core-collapse supernovae Woosley
and Bloom (2006) and short GRBs with neutron star mergers Nakar (2007). There
are however many exceptions to this conventional dichotomy, and many claims
for additional GRB sub-classes have also been discussed in literature Kulkarni
and Desai (2017); Bhave et al. (2022) (and references therein).

GRBs have also been proposed as distance indicators or standard candles over
the past two decades due to a large number of observed correlations between
various GRB observables in both the prompt and afterglow emission phase
Dainotti and Amati (2018); Parsotan and Ito (2022); Dainotti and Del Vecchio
(2017); Luongo and Muccino (2023); Dainotti et al. (2022); Pradyumna and Desai
(2022). These correlations have often been used to estimate cosmological
parameters Moresco et al. (2022). However, there is an inherent circularity
problem in using the GRB observables as cosmological probes, since the
calculation of some of these observables is based on an underlying
cosmological model. To get around this circularity problem, two methodologies
have been used in literature. One way is to simultaneously constrain the GRB
correlations along with the cosmological parameters Khadka et al. (2021).
Alternately, a large number of ancillary probes have also been used to obtain
cosmology-independent estimates of distances corresponding to the GRB redshift
such as Type Ia SN, cosmic chronometers, BAO $H(z)$ measurements, X-ray and UV
luminosities of quasars, galaxy clusters etc Govindaraj and Desai (2022) (and
references therein).

Here, we focus on correlations between six pairs of GRB observables:
$\tau_{lag}-L$, $V-L$, $E_{peak}-L$, $E_{peak}-E_{\gamma}$, $\tau_{RT}-L$,
$E_{peak}-E_{iso}$, which were first proposed in Wang et al. (2011). We note
that among the above correlations, the $E_{peak}-E_{\gamma}$ relation is often
refered to as as the Amati relation in literature Amati (2006). The
aforementioned work simultaneously fitted for cosmology and as well any
possible correlations Wang et al. (2011). These same correlations were then
analyzed by obtaining model-independent distances to GRBs using the Pantheon
compilation Scolnic et al. (2018) of Type Ia SN in Tang et al. (2021) (T21,
hereafter). Here, we carry out a variant of the analysis done in Tang et al.
(2021), by using X-ray and UV luminosities of quasars instead of Type Ia SN to
probe the same correlations first considered in Wang et al. (2011). We note
that quasars have previously being used to probe the efficacy of the Amati
relation in GRBs Dai et al. (2021) using the datasets in Khadka et al. (2021);
Demianski et al. (2017a, b).

The outline of this manuscript is as follows. We discuss the datasets used in
this work in Sect. II. Our analysis and results can be found in Sect. III. We
conclude in Sect. IV.

##  II Datasets

We use the same GRB datasets as those considered in Wang et al. (2011). This
work considered a sample of 116 long GRBs with redshifts between 0.17 to 8.2.
This sample consisted of GRBs observed by SWIFT until 2012 in conjunction with
additional GRBs from the pre-SWIFT era Schaefer (2007). The observables
assembled for each GRB consists of bolometric peak flux ($P_{bolo}$),
bolometric fluence ($S_{bolo}$), beaming factor ($F_{beam}$), time lag between
low and high energy photon light curves ($\tau_{lag}$), peak energy of the
spectrum ($E_{peak}$), minimum rise time of the peaks for which the light
curve rise by half its peak flux ($\tau_{RT}$), and the variability of the
light curve ($V$). More details can also be found in Wang et al. (2011). In
addition to the values, the $1\sigma$ error bars were also provided for each
of the above observables in Wang et al. (2011). Among these variables,
$\tau_{RT}$, $E_{peak}$, $\tau_{lag}$, and $V$ were obtained directly from the
light curves, whereas $P_{bolo}$ and $S_{bolo}$ can be obtained from the
observed GRB spectrum as outlined in Wang et al. (2011). In order to test for
potential correlations, additional quantities such as the isotropic peak
luminosity ($L$), isotropic equivalent energy ($E_{iso}$), collimation-
corrected energy ($E_{\gamma}$) depend on the luminosity distance ($d_{L}$),
which in turn need an underlying cosmological model. We now discuss the
quantities which depend on $d_{L}$. The relation between $L$ and $P_{bolo}$ is
given by:

$L=4\pi d_{L}^{2}P_{bolo},$ (1)

$E_{iso}$ is related to $D_{L}$ using:

$E_{iso}=4\pi d_{L}^{2}S_{bolo}(1+z)^{-1}$ (2)

Finally $E_{\gamma}$ is given by:

$E_{\gamma}=4\pi d_{L}^{2}S_{bolo}F_{beam},$ (3)

where $F_{beam}$ is the beaming factor which was estimated using the empirical
formula derived in Sari et al. (1999).

This dataset was used to study six different pairs of luminosity correlations:
$\tau_{lag}-L$, $V-L$, $E_{peak}-L$, $E_{peak}-E_{\gamma}$, $\tau_{RT}-L$,
$E_{peak}-E_{iso}$, where some of the above variables were scaled according to
redshift as explained in the next section. To circumvent this circularity
problem in calculating some of the above variables in a model agnostic manner,
a simultaneous fit to a linear regression between the above variables and the
underlying cosmology model was simultaneously done Wang et al. (2011). The
intrinsic scatter of the $V-L$ correlation was found to be very large, but the
other variables had a tight correlation with a negligible redshift evolution.
These same set of correlations between the variables were then re-examined
using a theory-agnostic approach without assuming an underlying cosmological
model Tang et al. (2021). For this purpose, a model-independent estimate of
$d_{L}$ at each GRB redshift was obtained using deep learning and Gaussian
process based regression using Type Ia supernova data. This work also tested
for a redshift evolution for the same six regression relations (considered in
Wang et al. (2011)) by dividing the GRB dataset into a low redshift and high
redshift sample. Among these, only $Ep-E_{\gamma}$ was found to contain no
redshift evolution Tang et al. (2021).

###  II.1 Quasar dataset

In order to obtain the distance corresponding to a given GRB redshift, instead
of Type Ia SN, we use quasars as distance anchors. Quasars contain Active
galactic nuclei, where the energy release occurs due to the accretion onto a
supermassive black hole. Quasars have been detected up to a redshift of $z=7$
and are therefore comparable to those found for GRBs. Quasars have been
observed throughout the electromagnetic spectrum Mortlock et al. (2011). A
tight scaling relation between the optical-UV flux at rest frame wavelength of
2500 Åand X-ray flux at rest frame energy of 2 keV of a quasar has been
asserted (with a scatter of around 0.12 dex), which is independent of redshift
Risaliti and Lusso (2019).

The luminosity distance ($D_{L}$) at a given redshift is obtained from the
quasar X-ray and UV flux as follows Risaliti and Lusso (2019); Lusso et al.
(2020):

$\log D_{L}=\frac{[\log F_{X}-\beta-\gamma(\log
F_{UV}+27.5)]}{2(\gamma-1)}-0.5\log(4\pi)+28.5,$ (4)

where $F_{X}$ and $F_{UV}$ are the flux densities (in erg/s/$cm^{2}$/Hz)
$\beta\approx 8.2$ and $\gamma=0.6$. This relation was obtained assuming a
Hubble constant value of $H_{0}=70$ km/sec/Mpc. Lusso et al Lusso et al.
(2020) constructed a clean sample of 2,421 optically selected quasars spanning
the redshift range $0.006\leq z\leq 7.52$ where distance modulus (DM) and
associated errors were obtained using $D_{L}$ from the above equation as
below:

$\rm{DM}=5\log(D_{L})-5\log(10pc)$ (5)

Therefore one can obtain the distance modulus for each quasar redshift from
the quasar X-ray and UV fluxes.

##  III Analysis and Results

The first step in the analysis involves obtaining a model-independent estimate
of $D_{L}$ from the quasar Dataset. As mentioned in the previous section, we
start with the distance modulus ($\mu$) and redshift ($z$) data for 2,421
quasars. With this data, we carry out a non-parametric reconstruction of
distance modulus ($\mu$) at any redshift $z$ using Gaussian Process Regression
(GPR). GPR is a generalization of a Gaussian distribution, characterized by a
mean and a covariance function (usually called the kernel function) Seikel et
al. (2012). More details about GPR can be found in our previous works
Singirikonda and Desai (2020); Bora and Desai (2021a, b); Bora et al. (2022);
Agrawal et al. (2021); Mendonça et al. (2021a, b); Holanda et al. (2022);
Govindaraj and Desai (2022). For the GPR implementation, we use the publicly
available Python package sklearn Pedregosa et al. (2011) to reconstruct the
distance modules ($\mu$) as a function of redshift ($z$) as follows:

k(x, x’) = ConstantKernel() + 1.0 * DotProduct(1) ** 0.1 (6)

The kernel we have used is a sum of linear and constant kernels. A linear
kernel with an exponent captures relations in the data, and a constant kernel
is used to scale magnitude. The GPR reconstruction of distance modules ($\mu$)
and the associated $1\sigma$ error bars can be found in Fig. 1. Once we have
reconstructed distance modulus ($\mu$) at any $z$, we can estimate the
luminosity distance ($D_{L}$) from $\mu$, by inverting Eq. 5. The errors in
$D_{L}$ can be obtained using standard error propagation. Once we obtain
$D_{L}$ for every GRB redshift, we can calculate $L$ and $E_{iso}$ and using
Eq. 1 and Eq. 2, respectively.

Figure 1: Reconstruction of the distance modulus of $\mu$ with GPR using X-ray
and UV fluxes for 2,421 quasars.

We then construct the following linear regression relations between the six
aforementioned GRB quantities in log space as follows:

$\displaystyle\log\frac{L}{\operatorname{erg}\mathrm{s}^{-1}}$
$\displaystyle=a_{1}+b_{1}\log\frac{\tau_{\operatorname{lag},i}}{0.1s},$ (7)
$\displaystyle\log\frac{L}{\operatorname{erg~{}s}^{-1}}$
$\displaystyle=a_{2}+b_{2}\log\frac{V_{i}}{0.02},$ (8)
$\displaystyle\log\frac{L}{\operatorname{erg~{}s}^{-1}}$
$\displaystyle=a_{3}+b_{3}\log\frac{E_{p,i}}{300\mathrm{keV}}$ (9)
$\displaystyle\log\frac{E_{\gamma}}{\text{ erg }}$
$\displaystyle=a_{4}+b_{4}\log\frac{E_{p,i}}{300\mathrm{keV}},$ (10)
$\displaystyle\log\frac{L}{\text{ erg s }}$
$\displaystyle=a_{5}+b_{5}\log\frac{\tau_{\mathrm{RT},i}}{0.1\mathrm{~{}s}},$
(11) $\displaystyle\log\frac{E_{\text{iso }}}{\text{ erg }}$
$\displaystyle=a_{6}+b_{6}\log\frac{E_{p,i}}{300\mathrm{keV}}$ (12)

where $\tau_{\mathrm{RT},i}=\tau_{\mathrm{RT}}/(1+z)$,
$\tau_{\operatorname{lag},i}=\tau_{\operatorname{lag}}/(1+z)$, $V_{i}=V(1+z)$
and $E_{p},i=E(1+z)$.

To get the best-fit parameters for each of the above equations, we apply the
D’Agostini’s likelihood which incorporates errors in both the ordinate and
abscissa D’Agostini (2005):

$\mathcal{L}\left(\sigma_{\mathrm{int}},a,b\right)\propto\prod_{i}\frac{1}{\sqrt{\sigma_{\mathrm{int}}^{2}+\sigma_{yi}^{2}+b^{2}\sigma_{xi}^{2}}}\times\exp\left[-\frac{\left(y_{i}-a-bx_{i}\right)^{2}}{2\left(\sigma_{\mathrm{int}}^{2}+\sigma_{yi}^{2}+b^{2}\sigma_{xi}^{2}\right)}\right],$
(13)

where $x$ and $y$ denote the abscissa and ordinate and $\sigma_{xi}$ and
$\sigma_{yi}$ are the corresponding errors; $\sigma_{\mathrm{int}}$ denotes
the intrinsic scatter in each regression relation. To get the best-fit values
of each of the parameters, we do Bayesian regression and sample the posterior
using the emcee MCMC sampler Foreman-Mackey et al. (2013).

In order to test if the aforementioned correlations vary with redshift,
similar to Tang et al. (2021), we bifurcate the GRB datset into two subsamples
corresponding to the following redshift bins: the low-z sample ($z\leq 1.4$)
which consists of 50 GRBs, and the high-z sample ($z>1.4$) which consists of
66 GRBs. We investigate the redshift dependence of luminosity correlations for
these two subsamples. We also show the results for the full GRB sample.

The best-fit intervals at 68%, 90%, and 99% credible intervals for all six
regression relations can be found in Fig. 2 for the full GRB sample as well as
after bifurcating the sample based on redshift. The scatter plots between the
six pairs of variables along with the best-fits are shown in Fig. 3. We find
that our results for the slope, intercept and intrinsic scatter for almost all
the six scaling relations are consistent with that obtained in T21. The only
exception is the $E_{p}-E_{\gamma}$ or Amati correlation, where we see an
intrinsic scatter of 26-40%, whereas the intrinsic scatter was found to be
$<23\%$ in T21. We find that the values for the slope and intercept are
consistent or at most $1-2\sigma$ discrepant between the low redshift and high
redshift samples indicating negligible evolution of the scaling relations.
Furthermore, all the six regression relations show a high intrinsic scatter of
greater than 30%. This includes the Amati relation, where we get a large
scatter of 47%, when we consider the full data sample, similar to our results
of using galaxy clusters as distance anchors Govindaraj and Desai (2022).
Therefore, the regression relations between these observables cannot be used
as model-independent probes of cosmological parameters.

Figure 2: The 68%,90%, and 99% credible intervals along with the marginalized
PDFs for the parameters for the regression relations between six pairs of GRB
observables for low redshift ($z<1.4$), high redshift ($z>1.4$) as well as the
full GRB dataset.

Figure 3: The luminosity correlations for low-$z$ (blue dots) and high-$z$ (red dots) GRBs. Error bars denote the 1$\sigma$ uncertainties. The lines depict the best-fit results, blue line for low-$z$ GRBs, red line for high-$z$ GRBs and black line for all-$z$ GRBs. Table 1: The best-fitting parameters of GRB luminosity correlations. $N$ is the number of GRBs in each subsample. Correlation | Sample | N | $a$ | $b$ | $\sigma_{\rm int}$
---|---|---|---|---|---
$\tau_{\rm lag}-L$ | low-z | 37 | 51.97$\pm$0.09 | -0.77$\pm$0.14 | 0.44$\pm$0.07
| high-z | 32 | 52.14$\pm$0.07 | -0.59$\pm$0.12 | 0.35$\pm$0.06
| All-z | 69 | 52.04$\pm$0.06 | -0.69$\pm$0.09 | 0.39$\pm$0.04
$V-L$ | low-z | 47 | 52.00$\pm$0.24 | 0.63$\pm$0.36 | 0.88$\pm$0.13
| high-z | 57 | 52.41$\pm$0.18 | 0.25$\pm$0.17 | 0.63$\pm$0.07
| All-z | 104 | 52.19$\pm$0.12 | 0.43$\pm$0.13 | 0.73$\pm$0.06
$E_{p}-L$ | low-z | 50 | 51.75$\pm$0.09 | 1.38$\pm$0.18 | 0.58$\pm$0.07
| high-z | 66 | 51.99$\pm$0.06 | 1.08$\pm$0.15 | 0.39$\pm$0.04
| All-z | 116 | 51.88$\pm$0.05 | 1.29$\pm$0.11 | 0.48$\pm$0.04
$E_{p}-E_{\gamma}$ | low-z | 12 | 50.55$\pm$0.09 | 1.39$\pm$0.20 | 0.26$\pm$0.09
| high-z | 12 | 50.30$\pm$0.14 | 1.37$\pm$0.44 | 0.39$\pm$0.16
| All-z | 24 | 50.43$\pm$0.07 | 1.39$\pm$0.17 | 0.29$\pm$0.07
$\tau_{\rm RT}-L$ | low-z | 39 | 52.52$\pm$0.12 | -1.24$\pm$0.18 | 0.46$\pm$0.06
| high-z | 40 | 52.38$\pm$0.08 | -0.75$\pm$0.17 | 0.42$\pm$0.06
| All-z | 79 | 52.42$\pm$0.07 | -1.02$\pm$0.11 | 0.43$\pm$0.041
$E_{p}-E_{\rm iso}$ | low-z | 40 | 52.43$\pm$0.10 | 1.49$\pm$0.20 | 0.59$\pm$0.08
| high-z | 61 | 52.51$\pm$0.06 | 1.22$\pm$0.14 | 0.39$\pm$0.04
| All-z | 101 | 52.47$\pm$0.05 | 1.38$\pm$0.12 | 0.48$\pm$0.04

##  IV Conclusions

In a recent work, T21 tested the empirical correlations among six pairs of GRB
observables for 116 long GRBs (first considered in Wang et al. (2011)) in
order to test their efficacy as a cosmological probe. However, one needs an
estimate of the luminosity distance, which depends on an underlying
cosmological model to calculate some of these GRB observables. In order to get
around this circularity problem, luminosity distances from Type Ia supernovae
were used as distance anchors, and the corresponding distance at a given GRB
redshift was obtained using ANN based regression.

In this work, we carry out the same exercise as T21, but use X-ray and UV
fluxes of quasars instead of Type Ia supernovae in order to get the luminosity
distance at a given GRB redshift. The interpolation has been done using
Gaussian process regression. Similar to Tang et al. (2021), we test the
correlations for both the low redshift and high redshift sample, after
bifurcating the dataset at $z=1.4$. Our results for the best-fit values for
all the six regression relations can be found in Table 1. The marginalized
credible intervals are shown in Fig. 2, whereas the scatter plots for the six
regression relations are shown in Fig. 3. Our conclusions are as follows:

  * •

Our slopes and intercepts agree with the corresponding results with T21 for
both the low and high redshift as well as the full sample.

  * •

Our intrinsic scatter for almost all the scaling relations are comparable to
that found in T21. The only exception is the Amati relation where we see a
much higher intrinsic scatter compared to T21.

  * •

Although there is negligible redshift evolution in the scaling relations, the
high intrinsic scatter implies that we cannot use these observables for model-
independent estimate of cosmological parameters.

## Acknowledgements

This work builds upon the M.Tech thesis of Shreeprasad Bhat and we are
grateful to him for sharing his codes.

## References

  * Tang et al. (2021) L. Tang, X. Li, H.-N. Lin, and L. Liu, Astrophys. J. 907, 121 (2021), eprint 2011.14040.
  * Kumar and Zhang (2015) P. Kumar and B. Zhang, Phys. Rep. 561, 1 (2015), eprint 1410.0679.
  * Yu et al. (2022) Y.-W. Yu, H. Gao, F.-Y. Wang, and B.-B. Zhang, in _Handbook of X-ray and Gamma-ray Astrophysics. Edited by Cosimo Bambi and Andrea Santangelo_ (2022), p. 31.
  * Singh and Desai (2022) A. Singh and S. Desai, JCAP 2022, 010 (2022), eprint 2108.00395.
  * Desai (2023) S. Desai, arXiv e-prints arXiv:2303.10643 (2023), eprint 2303.10643.
  * Kouveliotou et al. (1993) C. Kouveliotou, C. A. Meegan, G. J. Fishman, N. P. Bhat, M. S. Briggs, T. M. Koshut, W. S. Paciesas, and G. N. Pendleton, Astrophys. J. Lett. 413, L101 (1993).
  * Woosley and Bloom (2006) S. E. Woosley and J. S. Bloom, Ann. Rev. Astron. Astrophys. 44, 507 (2006), eprint astro-ph/0609142.
  * Nakar (2007) E. Nakar, Phys. Rep. 442, 166 (2007), eprint astro-ph/0701748.
  * Kulkarni and Desai (2017) S. Kulkarni and S. Desai, Astrophysics & Space Sciences 362, 70 (2017), eprint 1612.08235.
  * Bhave et al. (2022) A. Bhave, S. Kulkarni, S. Desai, and P. K. Srijith, Astrophysics & Space Sciences 367, 39 (2022), eprint 1708.05668.
  * Dainotti and Amati (2018) M. G. Dainotti and L. Amati, PASP 130, 051001 (2018), eprint 1704.00844.
  * Parsotan and Ito (2022) T. Parsotan and H. Ito, Universe 8, 310 (2022), eprint 2204.09729.
  * Dainotti and Del Vecchio (2017) M. G. Dainotti and R. Del Vecchio, New Astronomy Reviews 77, 23 (2017), eprint 1703.06876.
  * Luongo and Muccino (2023) O. Luongo and M. Muccino, Mon. Not. R. Astron. Soc. 518, 2247 (2023), eprint 2207.00440.
  * Dainotti et al. (2022) M. G. Dainotti, V. Nielson, G. Sarracino, E. Rinaldi, S. Nagataki, S. Capozziello, O. Y. Gnedin, and G. Bargiacchi, Mon. Not. R. Astron. Soc. 514, 1828 (2022), eprint 2203.15538.
  * Pradyumna and Desai (2022) S. Pradyumna and S. Desai, Journal of High Energy Astrophysics 35, 77 (2022), eprint 2204.03363.
  * Moresco et al. (2022) M. Moresco, L. Amati, L. Amendola, S. Birrer, J. P. Blakeslee, M. Cantiello, A. Cimatti, J. Darling, M. Della Valle, M. Fishbach, et al., Living Reviews in Relativity 25, 6 (2022), eprint 2201.07241.
  * Khadka et al. (2021) N. Khadka, O. Luongo, M. Muccino, and B. Ratra, JCAP 2021, 042 (2021), eprint 2105.12692.
  * Govindaraj and Desai (2022) G. Govindaraj and S. Desai, JCAP 2022, 069 (2022), eprint 2208.00895.
  * Wang et al. (2011) F.-Y. Wang, S. Qi, and Z.-G. Dai, Mon. Not. R. Astron. Soc. 415, 3423 (2011), eprint 1105.0046.
  * Amati (2006) L. Amati, Mon. Not. R. Astron. Soc. 372, 233 (2006), eprint astro-ph/0601553.
  * Scolnic et al. (2018) D. M. Scolnic, D. O. Jones, A. Rest, Y. C. Pan, R. Chornock, R. J. Foley, M. E. Huber, R. Kessler, G. Narayan, A. G. Riess, et al., Astrophys. J. 859, 101 (2018), eprint 1710.00845.
  * Dai et al. (2021) Y. Dai, X.-G. Zheng, Z.-X. Li, H. Gao, and Z.-H. Zhu, Astron. & Astrophys. 651, L8 (2021), eprint 2111.05544.
  * Demianski et al. (2017a) M. Demianski, E. Piedipalumbo, D. Sawant, and L. Amati, Astron. & Astrophys. 598, A112 (2017a), eprint 1610.00854.
  * Demianski et al. (2017b) M. Demianski, E. Piedipalumbo, D. Sawant, and L. Amati, Astron. & Astrophys. 598, A113 (2017b), eprint 1609.09631.
  * Schaefer (2007) B. E. Schaefer, Astrophys. J. 660, 16 (2007), eprint astro-ph/0612285.
  * Sari et al. (1999) R. Sari, T. Piran, and J. P. Halpern, Astrophys. J. Lett. 519, L17 (1999), eprint astro-ph/9903339.
  * Mortlock et al. (2011) D. J. Mortlock, S. J. Warren, B. P. Venemans, M. Patel, P. C. Hewett, R. G. McMahon, C. Simpson, T. Theuns, E. A. Gonzáles-Solares, A. Adamson, et al., Nature (London) 474, 616 (2011), eprint 1106.6088.
  * Risaliti and Lusso (2019) G. Risaliti and E. Lusso, Nature Astronomy 3, 272 (2019), eprint 1811.02590.
  * Lusso et al. (2020) E. Lusso, G. Risaliti, E. Nardini, G. Bargiacchi, M. Benetti, S. Bisogni, S. Capozziello, F. Civano, L. Eggleston, M. Elvis, et al., Astron. & Astrophys. 642, A150 (2020), eprint 2008.08586.
  * Seikel et al. (2012) M. Seikel, C. Clarkson, and M. Smith, JCAP 2012, 036 (2012), eprint 1204.2832.
  * Singirikonda and Desai (2020) H. Singirikonda and S. Desai, European Physical Journal C 80, 694 (2020), eprint 2003.00494.
  * Bora and Desai (2021a) K. Bora and S. Desai, JCAP 2021, 052 (2021a), eprint 2104.00974.
  * Bora and Desai (2021b) K. Bora and S. Desai, European Physical Journal C 81, 296 (2021b), eprint 2103.12695.
  * Bora et al. (2022) K. Bora, R. F. L. Holanda, S. Desai, and S. H. Pereira, European Physical Journal C 82, 17 (2022), eprint 2106.15805.
  * Agrawal et al. (2021) R. Agrawal, H. Singirikonda, and S. Desai, JCAP 2021, 029 (2021), eprint 2102.11248.
  * Mendonça et al. (2021a) I. E. C. R. Mendonça, K. Bora, R. F. L. Holanda, S. Desai, and S. H. Pereira, JCAP 2021, 034 (2021a), eprint 2109.14512.
  * Mendonça et al. (2021b) I. E. C. R. Mendonça, K. Bora, R. F. L. Holanda, and S. Desai, JCAP 2021, 084 (2021b), eprint 2107.14169.
  * Holanda et al. (2022) R. F. L. Holanda, K. Bora, and S. Desai, European Physical Journal C 82, 526 (2022), eprint 2105.10988.
  * Pedregosa et al. (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, et al., Journal of Machine Learning Research 12, 2825 (2011).
  * D’Agostini (2005) G. D’Agostini, arXiv e-prints physics/0511182 (2005), eprint physics/0511182.
  * Foreman-Mackey et al. (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publ. Astron. Soc. Pac. 125, 306 (2013), eprint 1202.3665.