#  Observational constraints on the EDSFD parametrization in $f(R,T)$ theory
of gravity

Harshna Balhara<EMAIL_ADDRESS>Department of Mathematics, Netaji
Subhas University of Technology, New Delhi-110 078, India  J. K. Singh
<EMAIL_ADDRESS>Department of Mathematics, Netaji Subhas University of
Technology, New Delhi-110 078, India  J. Jena<EMAIL_ADDRESS>Department of Mathematics, Netaji Subhas University of Technology, New
Delhi-110 078, India

###### Abstract

This paper investigates a cosmological model using a suitable parameterization
represented as a differential equation with respect to the energy density of
the scalar field, $\rho_{\phi}$, in the context of $f(R,T)$ gravity theory for
the flat Friedmann-Lema$\hat{\imath}$tre-Robertson-Walker (FLRW) model. This
parameterization is known as the Energy Density Scalar Field Differential
Equation (EDSFD) parametrization. It results in a solution of the Hubble
parameter containing four model parameters namely, $\Omega_{m0},\Omega_{\phi
0},H_{0},$ and $\alpha$. To constrain the model parameters, $77$ data points
from the Hubble dataset, $1048$ points from the Pantheon dataset, and BAO of
$6$ data points are used. We analyze and compare our model with the standard
$\Lambda$CDM model using the constrained values. The evolution of the physical
parameters containing the deceleration parameter, density parameter, Equation
of State (EoS) for Dark Energy, and $Om(z)$ diagnostic has been discussed. It
has been observed that our model exhibits quintessence-like behavior in the
redshift range $-1<z<1$ and tends to $\Lambda$CDM in the infinite future.

PACS numbers: 04.50.-h, 98.80.-k.

Keywords: Modified theory, EDSFD parametrization, Observational constraints.

##  I  Introduction

It is now generally accepted that the expansion of the cosmos at redshifts
lower than unity underwent a ”dynamic phase change” from a decelerating to an
accelerating state, which has been confirmed by numerous independent analyses
[1, 2, 3]. The measurements of the luminosity distance versus redshift
relation using supernovae (SNeIa) of type Ia, which were initially interpreted
in the context of $\Lambda$CDM scenarios using either background or
inhomogeneous luminosity distances, provide the most direct evidence for the
universe’s current accelerating stage [4, 5, 6, 7, 8]. Independent theoretical
observational analyses, however, indicate more generic models with a negative-
pressure dark energy component as their main component [9, 10, 11].

The current cosmic picture strongly supports the following: the spatial
geometry (curvature) is flat, and the dynamics are controlled by a component
known as dark energy (DE), which makes up $3/4^{th}$ of the composition and
$1/4^{th}$ of the matter component (baryons plus dark). The simplest and
theoretically most appealing way to characterize this dark energy component is
by using a positive cosmological constant $\Lambda$. The cosmological constant
$\Lambda$ has an equation of state (EoS), $\omega=\frac{p}{\rho}=-1$, where
$p$ and $\rho$ are the effective thermodynamic pressure and energy density
associated with $\Lambda$. The dynamical DE scenario is an alternate DE
proposal. Numerous scalar field dark energy models have been discussed so far,
including the quintessence [12, 13, 14, 15], K-essence [16], tachyon [17],
Phantom [18], ghost condensate [19, 20], quintom [21], braneworld models [22],
and Chaplygin gas models [23] etc.

Modified theories of gravity have received a lot of attention lately,
especially in light of their potential contribution to the universe’s
accelerating expansion. As a result of its simplicity and lack of ghosts, the
so-called $f(R)$ gravity [24, 25, 26, 27, 28, 29] is the simplest and most
widely accepted classical extension of General Relativity (GR). In addition,
compared to previous modifications of gravity, the reconstruction of $f(R)$
theories that can replicate the dark energy epoch and even the inflationary
phase is typically simple [30]. It is possible for $f(R)$ gravity to replicate
the $\Lambda$CDM era, imitate a cosmological constant in the present, and even
unify all of cosmological history. Additionally, some $f(R)$ gravities, or so-
called viable models, can replicate a plausible cosmic evolution and pass the
local gravitational criteria [31]. Studying $f(R)$ gravity models allows us to
explore the connection between early-time inflation and late-time acceleration
[32, 33]. Existing literature highlights the crucial role played by higher-
order curvature terms in $f(R)$ gravity models in preventing cosmological
singularities [34, 35, 36].

In their work [37], Nojiri and Odinstov examined a modified gravity theory
that incorporates terms proportional to $ln(R)$ or $R^{-n}(lnR)^{m}$, which
exhibit growth at low curvatures. These $ln(R)$ or $R^{-n}(lnR)^{m}$ terms in
$f(R)$ gravity models are believed to contribute to the acceleration of the
universe. Nojiri and Odinstov further investigated the cosmology of $f(R)$
gravity, specifically considering the form $f(R)=R+\gamma
R^{-n}ln(\frac{R}{\mu^{2}})^{m}$ [38]. These functional forms were also
utilized in previous studies [37, 38, 39] to explore various aspects of the
theory. Additionally, in [40], it was demonstrated that all these models
exhibit the current accelerating phase of the universe, and the duration of
this phase relies on the coupling constants of the gravitational action. The
models of the type $f(R)=R-\frac{\alpha}{R^{n}}$, encounter challenges in
successfully passing the solar system tests [41] and are susceptible to
gravitational instabilities [42]. Furthermore, these theories fail to generate
a standard matter-dominated era preceding accelerated expansion [43, 44]. On
the other hand, the $f(R)=R+\frac{\alpha}{R^{n}}+\frac{\beta}{R^{m}}$ type
models face difficulties in meeting the range of constraints imposed by early
and late-time acceleration, big bang nucleosynthesis, and fifth-force
experiments [45].

Considering the numerous challenges associated with the majority of $f(R)$
models mentioned above, more complicated theories, particularly those that
have an unconventional coupling of curvature with the energy-momentum tensor
were explored. A family of modified gravity theories that contain the general
function $f(R,T)$ in the gravitational action, where $R$ and $T$ stand for the
curvature and trace of the energy-momentum tensor, respectively, has recently
been presented [46]. Since the inception of this theory, extensive research
has been conducted on various facets, including the thermodynamic properties
[47, 48, 49], energy conditions [50, 51], cosmological solutions derived from
analyzing a homogeneous and isotropic spacetime using phase-space analysis
[52], anisotropic cosmology [53, 54], cosmological solutions involving an
auxiliary scalar field [55] and investigation of additional aspects [56, 57,
58]. Studies on the cosmological evolution of $f(R,T)$ theory have taken into
account future singularities and the reconstruction of cosmological solutions
[59, 60, 61]. However, such a tight coupling between the curvature and the
trace $T$ implies a violation of the common continuity equation, a problem
that can be resolved by a suitable function $f(R,T)$ as first demonstrated in
Ref [62]. It is commonly understood that theories of modified gravity can be
classified into two distinct categories depending on the relationship between
matter and geometry, namely minimal coupling and non-minimal coupling. There
are several studies that examine the cosmological consequences of the non-
minimal coupling of matter and geometry in $f(R,T)$ gravity. Shabani and
colleagues [52] conducted a dynamical systems analysis on various Lagrangians
in modified gravity, namely a minimal $g(R)+h(T)$ Lagrangian, a pure non-
minimal $g(R)h(T)$ Lagrangian, and a non-minimal $g(R)(1+h(T))$ Lagrangian,
using the FLRW metric as a background. Several studies explore the
cosmological consequences of the non-minimal coupling of matter and geometry
in $f(R,T)$ gravity [63, 64, 51]. Numerous studies can be found in the
literature that focus on specific linear versions of $f(R,T)$ gravity [50, 65,
66, 67, 68]. The focus of our study is on $f(R,T)$ gravity with minimal
coupling, specifically $f(R,T)=R+2\lambda T$, where $\lambda$ represents an
arbitrary coupling constant.

The following is the presentation and structure of the paper: Section II
provides a brief overview of the formulation of field equations and its
solution in $f(R,T)$ gravity. In this section, by employing an appropriate
parameterization expressed as a differential equation involving the energy
density of the scalar field, $\rho_{\phi}$, we derive equations that
characterize the Hubble and deceleration parameters in terms of the redshift.
Section III provides a concise overview of the observational Hubble, Pantheon,
and BAO data that are utilized to assess the cosmological model and establish
limits on the model parameters. In this context, our aim is to demonstrate the
close correspondence between observational and theoretical outcomes. In
Section IV, we deliberate the outcomes obtained from the observational data
and observe the evolution of the universe. This section delves into the
examination of the jerk and snap parameters, revealing that our current model
aligns with a dark energy quintessence model. Additionally, we have depicted
the behavior of the densities and pressure along with the Equation of State
(EoS) parameter. The viability of our model is elucidated by our explanation
of energy conditions. Finally, we study the thermodynamical behavior of our
model in Section V and conclude our work in Section VI.

##  II Field Equations

Among the modified theories that are based on nonminimal interaction between
matter and geometry, the $f(R,T)$ gravity is an intriguing contender. This
updated theory’s behavior takes the form

$S=\frac{1}{\kappa^{2}}\int[f(R,T)+\mathcal{L}]\sqrt{-g}dx^{4}.$ (1)

Here $\kappa^{2}=1$ is the gravitational coupling constant, $R$ is the Ricci
scalar, $T$ is the trace of the stress energy-momentum tensor, while
$\mathcal{L}=\mathcal{L}_{m}+\mathcal{L}_{\phi}$ where $\mathcal{L}_{m}$ is
the matter Lagrangian and $\mathcal{L}_{\phi}$ is the action for the scalar
field. The action of the scalar field is given by

$\mathcal{L}_{\phi}={\Big{[}-\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi-V(\phi)\Big{]}\sqrt{-g}d^{4}x}.$
(2)

We analyze the evolution and dynamics of the universe for a particular form of
$f(R,T)$, namely $f=f_{1}(R)+f_{2}(T)$ where $f_{1}(R)=R$ and
$f_{2}(T)=2\lambda T$, $\lambda$ being an arbitrary coupling constant. We
assume that the universe is isotropic and homogeneous, with the energy
density, $\rho$, the thermodynamic pressure, $p$, and the matter Lagrangian
$\mathcal{L}_{m}=-\rho_{m}$. The Friedmann-Lema$\hat{\imath}$tre-Robertson-
Walker (FLRW) formula which explains the geometry of space-time is given by

$ds^{2}=dt^{2}-a^{2}(t)(dx^{2}+dy^{2}+dz^{2}),$ (3)

where $a(t)$ is the scale factor of the universe. The gravitational field
equations in $f(R,T)$ gravity can be found by varying the action, given by (1)
with respect to $g_{\mu\nu}$ as

$G_{\mu\nu}+(g_{\mu\nu}\square-\nabla_{\mu}\nabla_{\nu})=[8\pi+2f^{\prime}(T)]T_{\mu\nu}+2[f^{\prime}(T)p+\frac{1}{2}f(T)]g_{\mu\nu},$
(4)

The energy-momentum tensor is given by [69]

$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}\mathcal{L})}{\delta
g^{\mu\nu}},$ (5)

By assuming a matter action that solely depends on the metric and not on its
first derivatives, the energy-momentum tensor thus produces,

$T_{\mu\nu}=g_{\mu\nu}\mathcal{L}-2\frac{\delta\mathcal{L}}{\delta
g^{\mu\nu}}.$ (6)

The independent cosmological field equations for the FLRW metric are

$3H^{2}=(8\pi+3\lambda)\rho_{total}-\lambda p_{total},$ (7)
$2\dot{H}+3H^{2}=-(8\pi+3\lambda)p_{total}+\lambda\rho_{total}.$ (8)

where we take $8\pi=1$, $\rho_{total}=\rho_{m}+\rho_{\phi}$,
$p_{total}=p_{\phi}$ and dot indicates differentiation w.r.t. time $t$.
$\rho_{m}$, $\rho_{\phi}$, and $p_{\phi}$ denote the energy density of matter,
energy density of scalar field, and pressure of scalar field respectively. The
isotropy and homogeneity of the model can be explained by the dominant
component of the scalar field $\phi$, which only varies with time. Since the
scalar field $\phi$, is time-dependent, it can be treated as a perfect fluid
with energy density, $\rho_{\phi}$, and pressure, $p_{\phi}$. Assuming that
the scalar field, $\phi$, is the sole source of DE, with a potential $V(\phi)$
that interacts with itself, we can consider energy density, $\rho_{\phi}$, and
pressure, $p_{\phi}$, as the canonical components of the scalar field $\phi$,
within the framework of FLRW cosmology.

$\rho_{\phi}=\frac{\dot{\phi}^{2}}{2}+V(\phi),\qquad
p_{\phi}=\frac{\dot{\phi}^{2}}{2}-V(\phi).$ (9)

In general relativity, the conservation of the energy-momentum tensor is
expressed as follows [70]:

$\dot{\rho}_{m}+3\bigg{(}\frac{\dot{a}}{a}\bigg{)}\rho_{m}=0,$ (10)
$\dot{\rho}_{\phi}+3\bigg{(}\frac{\dot{a}}{a}\bigg{)}(\rho_{\phi}+p_{\phi})=0.$
(11)

Eq. (10) yields

$\rho_{m}=\rho_{m_{0}}a^{-3},$ (12)

where $\rho_{m_{0}}$ is an arbitrary integration constant. Solving Eq. (11),
we get

$\dot{\rho}_{\phi}=-3\frac{\dot{a}}{a}(1+\omega_{\phi})\rho_{\phi}.$ (13)

The EoS parameter of the scalar field $\phi$ is expressed as
$\omega_{\phi}=\frac{p_{\phi}}{\rho_{\phi}}$ and its value can be determined
using Eq. (11) as,

$\omega_{\phi}=-1-\frac{1}{3}a\frac{\rho_{\phi}^{\prime}}{\rho_{\phi}}.$ (14)

where prime denotes the derivative w.r.t. scale factor $a$. The system of Eqs.
(7), (8), (10) and (11) that only has three independent equations and four
unknowns $H$, $\rho_{m}$, $\rho_{\phi}$ and $p_{\phi}$ needs some additional
constraint equation to be solved. Here, we assume a suitable parametrization
in the form of a differential equation in $\rho_{\phi}$ as

$\rho_{\phi}^{\prime}+\alpha\Big{[}1-\frac{1}{\alpha(1+a)\log(1+a)}\Big{]}\rho_{\phi}=0,$
(15)

The generic answer to the differential Eq. (15) mentioned above is given by

$\rho_{\phi}(z)=e^{-\alpha a}\log(1+a),$ (16)

where $\alpha\in(0,1)$ is the model parameter. The energy density of scalar
field $\rho_{\phi}$ can be expressed in terms of redshift $z$ by using the
relationship between redshift $z$ and scale factor $a$, which is given by
$\frac{a}{a_{0}}=\frac{1}{1+z}$, where $a_{0}=1$ represents the current value
of the scale factor.

$\rho_{\phi}(z)=e^{\frac{-\alpha}{1+z}}\log\bigg{(}1+\frac{1}{1+z}\bigg{)},$
(17)

and

$\rho_{\phi_{0}}=e^{-\alpha}\log(2),$ (18)

Eqs. (17) and (18) yield

$\rho_{\phi}(z)=\frac{\rho_{\phi_{0}}}{\log(2)}e^{\frac{\alpha
z}{1+z}}\log\bigg{(}1+\frac{1}{1+z}\bigg{)}.$ (19)

where $\rho_{\phi_{0}}$ is the current value of the energy density of the
scalar field. Using Eq. (12), it is possible to express the energy density of
the matter field $\rho_{m}$ in terms of redshift $z$ as

$\rho_{m}(z)=\rho_{m_{0}}(1+z)^{3}.$ (20)

Combining Eqs. (7), (17) and (19) we obtain

$3H(z)^{2}=(3\lambda+1)\left((z+1)^{3}\rho_{\text{m0}}+\frac{\rho_{\text{$\phi$0}}e^{\frac{\alpha
z}{z+1}}\log\left(\frac{1}{z+1}+1\right)}{\log(2)}\right)\\\
-\lambda\left(-\frac{\frac{\rho_{\text{$\phi$0}}e^{\frac{\alpha
z}{z+1}}\left(\frac{\alpha}{z+1}-\frac{\alpha
z}{(z+1)^{2}}\right)\log\left(\frac{1}{z+1}+1\right)}{\log(2)}-\frac{\rho_{\text{$\phi$0}}e^{\frac{\alpha
z}{z+1}}}{(z+1)^{2}\left(\frac{1}{z+1}+1\right)\log(2)}}{3(z+1)}-\frac{\rho_{\text{$\phi$0}}e^{\frac{\alpha
z}{z+1}}\log\left(\frac{1}{z+1}+1\right)}{\log(2)}\right),$ (21)

The density parameter $\Omega=\rho/\rho_{c}$ is used to describe the contents
of the universe, with $\rho_{c}=3H^{2}/(8\pi G)^{2}$ representing the critical
density. In normalized units, $8\pi G=1$. Expressing Eq. (21) in terms of the
density parameter of matter and scalar field as

$H(z)^{2}=H_{0}^{2}(3\lambda+1)\left((z+1)^{3}\Omega_{\text{m0}}+\frac{\Omega_{\text{$\phi$0}}e^{\frac{\alpha
z}{z+1}}\log\left(\frac{1}{z+1}+1\right)}{\log(2)}\right)\\\
-H_{0}^{2}\lambda\left(-\frac{\frac{\Omega_{\text{$\phi$0}}e^{\frac{\alpha
z}{z+1}}\left(\frac{\alpha}{z+1}-\frac{\alpha
z}{(z+1)^{2}}\right)\log\left(\frac{1}{z+1}+1\right)}{\log(2)}-\frac{\Omega_{\text{$\phi$0}}e^{\frac{\alpha
z}{z+1}}}{(z+1)^{2}\left(\frac{1}{z+1}+1\right)\log(2)}}{3(z+1)}-\frac{\Omega_{\text{$\phi$0}}e^{\frac{\alpha
z}{z+1}}\log\left(\frac{1}{z+1}+1\right)}{\log(2)}\right),$ (22)

where $\Omega_{m_{0}}=\frac{\rho_{m_{0}}}{3H_{0}^{2}}$ and
$\Omega_{\phi_{0}}=\frac{\rho_{\phi_{0}}}{3H_{0}^{2}}$ are the present values
of matter and scalar field density parameters respectively.

$H=\frac{1}{(z+1)^{3/2}\sqrt{z+2}\sqrt{\log(8)}}\Bigg{(}\surd\Big{(}H_{0}^{2}(4.28365+z(27.8437+z(77.1057+z(117.8+z(107.091+z(57.8293+z\\\
(17.1346+2.14182z)))))))\Big{)}\Omega_{m_{0}}+\Omega_{\phi_{0}}(-0.01z-0.01)e^{\frac{\alpha
z}{z+1}}\\\ +e^{\frac{\alpha
z}{1+z}}(0.06+H_{0}^{2}(6.18+z(21.63+z(27.81+z(15.45+3.09z)))))\\\
+0.02\alpha+z(0.01\alpha+z((0.03z+0.15)z+0.27)+0.21))\Omega_{\phi_{0}}\log\left(\frac{1}{z+1}+1\right)\Bigg{)}.$
(23)

The following section involves constraining the model parameters
$\Omega_{m_{0}},\Omega_{\phi_{0}},H_{0}$, and $\alpha$ by utilizing
observational datasets and determining their optimal fit values. Subsequently,
these values will be utilized to analyze the behavior of various physical
parameters.

##  III Observation Dataset

In the previous section, we examined $f(R,T)$ gravity and found an exact
solution for the resulting field equations. The solution includes four model
parameters, namely, $\Omega_{m0},\Omega_{\phi 0},H_{0},$ and $\alpha$. In
order to verify our method, we need to limit these model parameters using
observational datasets to determine their best-fit values. Three datasets were
used in this study: the Hubble dataset, which includes 77 data points, the
recently released Pantheon sample, which includes 1048 Supernova Type Ia
experiment findings from surveys such as the Low-z, SDSS, SNLS, Pan-
STARRS1(PS1) Medium Deep Survey, and HST [71] in the redshift range
$z\in(0.01,2.26)$ and the BAO dataset consisting of 6 data points.

Figure 1:  The graph provides the best-fit values of the model parameters
$\Omega_{m0}$, $\Omega_{\phi 0}$, $H_{0},$ and $\alpha$ with $1-\sigma$ and
$2-\sigma$ errors for the 77-point Hubble dataset.

###  III.1 Hubble Data

Observational cosmology employs the Hubble parameter to examine the expansion
of the universe. Two widely used techniques for estimating $H(z)$ at certain
redshifts are the extraction of $H(z)$ from Baryon Acoustic Oscillation (BAO)
data and the differential age method. Ref. [72] provides a brief summary of
the revised dataset containing 77 points in the redshift range of
$0.07<z<2.42$. To calculate the model parameters, we utilize MCMC simulation
and the $\chi^{2}$ test. The value of $\chi^{2}_{Hubble}$ for the
observational Hubble parameter can be computed as

$\chi^{2}_{Hubble}=\sum_{i}^{77}\frac{[H_{th}(z_{i},\Omega_{m0},\Omega_{\phi
0},H_{0},\alpha)-H_{obs}(z_{i})]^{2}}{\sigma^{2}_{(z_{i})}},$ (24)

Here, the theoretical value of the Hubble parameter is denoted by $H_{th}$,
while $H_{obs}$ denotes the observed value. The standard error of the observed
Hubble parameter is represented by $\sigma_{(z_{i})}$.

Using the Hubble dataset consisting of $77$ points as mentioned earlier, we
generated two-dimensional contour plots with $1-\sigma$ and $2-\sigma$ errors
in Fig. 1 to determine the optimal values of the model parameters
$\Omega_{m0}$, $\Omega_{\phi 0}$, $H_{0},$ and $\alpha$. The best fit values
are obtained as $\Omega_{m0}=0.2542^{+0.0024}_{-0.0037}$, $\Omega_{\phi
0}=0.644^{+0.025}_{-0.029}$, $H_{0}=67.44^{+0.24}_{-0.34}$ and
$\alpha=0.184^{+0.053}_{-0.065}$.

###  III.2 Pantheon Data

The most recently published dataset for supernovae type $Ia$ is the Pantheon
sample, consisting of 1048 data points. In our analysis, we utilize these data
points, which have been confirmed spectroscopically by SNe$Ia$ and cover the
redshift range of $0.01<z<2.26$. The $\chi^{2}_{Pantheon}$ function for the
Pantheon dataset is taken to be,

$\chi_{Pantheon}^{2}=\sum\limits_{i=1}^{1048}\left[\frac{\mu_{th}(\mu_{0},\Omega_{m0},\Omega_{\phi
0},H_{0},\alpha)-\mu_{obs}(z_{i})}{\sigma_{\mu(z_{i})}}\right]^{2},$ (25)

Figure 2:  The plot depicts the best-fit values of the model parameters
$\Omega_{m0}$, $\Omega_{\phi 0}$, $H_{0},$ and $\alpha$ with $1-\sigma$ and
$2-\sigma$ errors for the $1048$ pantheon dataset.

In addition, the symbol $\sigma_{\mu(z_{i})}$ represents the standard error in
the actual value of $\mu$. The distance modulus $\mu_{th}$ is defined as
$\mu^{i}_{th}=\mu(D_{L})=m-M=5log_{10}D_{L}(z)+\mu_{0}$, where $m$ and $M$
represent the apparent and absolute magnitudes, respectively, and the nuisance
parameter $\mu_{0}$ is determined by $\mu_{0}=5log(H_{0}^{-1}/Mpc)+25$. To
obtain the luminosity distance $D_{L}$, one can utilize the formula
$D_{L}(z)=(1+z)H_{0}\int\frac{1}{H(z^{*})}dz^{*}$, where $H(z)$ is a series
that is limited to the tenth term and can be roughly integrated to determine
the luminosity distance. The most suitable values of the model parameters
$\Omega_{m0}$, $\Omega_{\phi 0}$, $H_{0}$ and $\alpha$ have been established
by examining two-dimensional contour plots with $1-\sigma$ and $2-\sigma$
errors, using the Pantheon dataset mentioned above. The resulting best-fit
values based on the $1048$ Pantheon dataset are
$\Omega_{m0}=0.327^{+0.019}_{-0.0097}$ , $\Omega_{\phi
0}=0.677^{+0.014}_{-0.014}$, $H_{0}=70.25^{+0.13}_{-0.13}$ and
$\alpha=0.195^{+0.055}_{-0.055}$ as shown in Fig. 2.

###  III.3 BAO Dataset

BAO represents a significant cosmological event that emerged during the
infancy of the universe. This BAO scale acts as a universal benchmark within
the vast structure of the cosmos, offering an unparalleled and resilient
cosmic measuring tool for gauging the historical expansion of the universe. In
this study, we have incorporated six points of the Baryon Acoustic Oscillation
(BAO) dataset, sourced from various surveys. These datasets encompass
observations from the Six Degree Field Galaxy Survey (6dFGS), the Sloan
Digital Sky Survey (SDSS), and the LOWZ samples of the Baryon Oscillation
Spectroscopic Survey (BOSS) [73, 74]. The expression for the scale of
dilation, denoted as $D_{v}(z)$, as presented in [75], is as follows:

$D_{v}(z)=\Big{(}\frac{d^{2}_{A}(z)z}{H(z)}\Big{)}^{1/3}$ (26)

here $d_{A}(z)$ represents the comoving angular diameter distance and is
formally defined as follows:

$d_{A}(z)=\int_{0}^{z}\frac{dz^{\prime}}{H(z^{\prime})},$ (27)

The chi-squared function for the BAO analysis is defined as

$\chi^{2}_{BAO}=A^{T}C^{-1}_{BAO}A,$ (28)

A’s dependency relies on the specific survey being considered, while $C^{-1}$
represents the inverse of the covariance matrix [76].

Figure 3:  The graph gives the best-fit values of the model parameters
$\Omega_{m0}$, $\Omega_{\phi 0}$, $H_{0},$ and $\alpha$ with $1-\sigma$ and
$2-\sigma$ errors for the 6 points of BAO dataset.

The optimal parameter values for the model, including $\Omega_{m0}$,
$\Omega_{\phi 0}$, $H_{0}$, and $\alpha$, have been determined by analyzing
two-dimensional contour plots, which incorporate $1-\sigma$ and $2-\sigma$
uncertainties, using the aforementioned BAO dataset. The resulting parameter
estimates based on the 6 data points from the BAO dataset are as follows:
$\Omega_{m0}=0.2105^{+0.0061}_{-0.0085}$, $\Omega{\phi
0}=0.607^{+0.078}_{-0.059}$, $H_{0}=70.44^{+0.12}_{-0.13}$, and
$\alpha=0.331^{+0.027}_{-0.033}$, as illustrated in Fig. 3.

Figure 4:  The contour plots for the model parameters $\Omega_{m0}$,
$\Omega_{\phi 0}$, $H_{0},$ and $\alpha$ with $1-\sigma$ and $2-\sigma$ errors
which provide the best-fit values of the model parameters from the joint
datasets of Hubble, Pantheon, and BAO datasets.

The left panel of Fig. 5 depicts the development of the Hubble function $H(z)$
over varying redshift $z$, where the obtained model is represented by colorful
lines. The blue bars on the curve display the 77 data points of the Hubble
dataset, along with their respective error bars. Additionally, the black
dotted line illustrates the $\Lambda$CDM model. The colored lines in the right
panel of Fig. 5 illustrate the relationship between $\mu(z)$ and $z$ for our
model, while the black dotted line represents the curve for the $\Lambda$CDM
model. Both models display a good match to the $H(z)$ data and Union 2.1
compilation data and their respective error bars.

(a)  (b)

Figure 5:  The error bar plots show the variations of the Hubble parameter
$H(z)$ and the distance modulus $\mu(z)$ in comparison to the $\Lambda$CDM.

Having obtained both theoretical formulas and numerical values for the model
parameters, we are now equipped to delve into the physical dynamics of the
model. Thus, the following section will be focused on exploring the physical
dynamics of additional significant cosmological parameters.

##  IV  Physical evaluation of the model

Numerous dark energy theories have been put forth in the literature to discuss
the accelerated expansion of the universe. Among various factors that trace
the dynamics of the universe, the deceleration parameter (DP), $q$, is said to
be one of the most important cosmic parameters [77]. The rate of acceleration
of the universe is described by $q$, which is a function of the Hubble
parameter and has the form $q=-1-\frac{\dot{H}}{H^{2}}$. If $q$ is positive,
the universe is decelerating, while a negative $q$ indicates an accelerated
phase. The values of $\Omega_{m0}$, $\Omega_{\phi 0}$, $H_{0},$ and $\alpha$
are utilized in the computation of $q$. Fig. 6a illustrates the behavior of
$q$ versus redshift $z$ and depicts the expansion from past to present. As is
evident, $q$ is in the accelerated stage at present which can be confirmed
from various cosmic observations that show the universe experiences a cosmic
acceleration in late periods with a slower pace of expansion [4, 7]. Negative
inflation in the universe is evidenced by the sudden halt in the expansion
after the phase transition during late times, approximately at $z=-1$ (refer
to Fig. 6a).

The Hubble, deceleration, and jerk parameters are purely kinematical, which
means they are not dependent on any gravity theory. These parameters are
solely related to the scale factor $a$ or redshift $z$ (as $a=1/1+z$). The
jerk parameter, which is a dimensionless third derivative of the scale factor
$a(t)$ w.r.t. cosmic time $t$, can be used as a straightforward way to look
for deviations from the standard $\Lambda$CDM model. It is defined as [78, 79,
80, 81, 82, 83]

$j=\frac{\dddot{a}}{aH^{3}}.$ (29)

The jerk parameter can be expressed in terms of the deceleration parameter as

$j=q(2q+1)+(1+z)\frac{dq}{dz}.$ (30)

In their study, Blandford et al. [78] explained how the jerk parameter can
serve as a convenient and alternate way of describing cosmological models that
are similar to the concordance $\Lambda$CDM model. The significant advantage
of using $j$ is that, for the $\Lambda$CDM model, $j$ always remains constant
and equal to $1$. It is worth mentioning that Sahni et al. [85, 84] emphasized
the significance of $j$ in distinguishing between various dark energy models,
as any deviation from $j=1$ would indicate a preference for a non-$\Lambda$CDM
model. Fig. 6b displays the model’s evolution of the jerk parameter, revealing
that for all numerically restricted values of the model parameters
$\Omega_{m0},\Omega_{\phi 0},H_{0},$ and $\alpha$, the value of $j_{0}$ for
our model falls within the range of 0.99 to 1.001. A value of approximately
$1$ for the cosmic jerk parameter at redshift $z=0$ indicates the occurrence
of cosmic acceleration.

(a)

(b)

Figure 6:  The behaviour of the deceleration parameter $q$ and jerk parameter
$j$ with respect to the redshift $z$.

In the study of cosmological models, the kinematic variables hold considerable
importance. Likewise, the equation of state (EoS) parameter, $\omega$, helps
in understanding the physical significance of energy sources during the
universe’s evolution. The EoS parameter is defined as

$\omega=\frac{p}{\rho}$ (31)

Based on observations, $\Lambda$CDM is a compelling contender for dark energy
with $\omega=-1$. Nonetheless, the universe’s acceleration stage is only
comprehensible when $\omega<-1/3$. This encompasses the quintessence phase,
where $0<\omega<-1$, while a value below $-1$ is referred to as the phantom
regime. The EoS parameter of the scalar field $\omega_{\phi}$ and the total
EoS parameter $\omega_{eff}$ have been examined in our model. The profile of
the EoS parameter of scalar field $\omega_{\phi}$ in relation to redshift $z$
is shown in the left panel of Fig. 7. It is noticeable that $\omega_{\phi}$
begins in the quintessence region during the universe’s evolution, remains in
that region up to the present time, and then experiences a sudden collapse
before converging to $\Lambda$CDM in the future. The quintessence phase can be
observed in the profile of $\omega_{eff}$ shown in the right panel of Fig. 7.
It is evident that $\omega_{eff}$ is less than zero, and the present value of
$\omega_{eff}$ lies between $-0.53$ and $-0.56$. This finding is consistent
with certain investigations [86], which suggest a period of acceleration at
present. In the future, it experiences a sudden collapse before emerging again
at later times.

(a)

(b)

Figure 7:  Evolution of the EoS parameter $\omega_{\phi}$ and $\omega_{eff}$
as a function of redshift $z$ for distinct values of model parameters
respectively.

We plot graphs illustrating the energy density as a function of the redshift
$z$, taking into account various constraint values of the model parameters
that were derived from different observational datasets. As the universe
evolves, the energy density is observed to be high initially but eventually
decreases such that $\rho_{eff}$ approaches zero in late times. This suggests
that the density of dark energy reduces to almost zero during later periods
(See Fig. 8b). The energy density of scalar field $\rho_{\phi}$ increases with
time but eventually collapses to $0$ at $z=-1$ (See Fig. 8a). The image in
Fig. 8c depicts isotropic pressure, $p_{eff}$, which is highly negative and
decreases with time. According to standard cosmology, cosmic acceleration is
associated with negative values of cosmic pressure.

Energy Conditions (ECs) are viewed as a sound generalization to the entire EMT
of the notion that the energy density of the universe can never be negative.
The use of ECs in examining the viability of numerous significant space-time
singularities problems, such as those involving black holes and wormholes, is
widespread. There are numerous ECs that are frequently employed in GR and
whose viabilities can be examined using the well-known Raychaudhuri equation.
The energy conditions can generally be expressed in one of two ways: (i)
geometrically, where the ECs are well expressed in terms of the Ricci tensor
or the Weyl tensor, or (ii) physically, where the ECs are expressed in terms
of the stress energy-momentum tensor or in terms of the energy density, which
functions as the time-like component and pressures $p_{i}$, i = 1, 2, 3, which
represent the three-space-like component. The typical point-wise ECs in GR are
defined as (a) Null Energy Condition : NEC: $\rho+p_{i}\geq 0,\forall i$;
Strong Energy Condition : SEC: $\rho+\sum_{i=1}^{3}p_{i}\geq 0,\forall i$;
Dominant Energy Condition : DEC: $\rho\geq 0,\rho\geq|p_{i}|,\forall i,$ where
$i=$ 1,2,3; Weak Energy Condition :WEC: $\rho\geq 0$, $\rho+p_{i}\geq
0,\forall i$.

(a)

(b)

(c)

Figure 8:  Evolution of $\rho_{\phi}$, $\rho_{eff}$, and $p_{eff}$ vs redshift
$z$ for distinct values of model parameters respectively..

By examining Fig. 9, it can be seen that our model demonstrates accelerating
expansion in the universe as NEC and DEC are satisfied, but SEC is violated.
The violation of SEC suggests that there may be exotic matter present in the
universe that produces anti-gravitational effects which are responsible for
the accelerated expansion of the universe (See Fig. 9c).

(a)

(b)

(c)

Figure 9:  Graphical behaviors of NEC, DEC, and SEC respectively.

$Om(z)$ stands for Om diagnostic [87, 88]. It is utilized to differentiate
between the basic $\Lambda$CDM model and several dark energy models. Om is
easier to reconstruct from observational data since it just uses the first
derivative of the scale factor via the Hubble parameter. $Om(z)$ is defined as

$Om(z)=\frac{(\frac{H(z)}{H_{0}})^{2}-1}{z(z^{2}+3z+3)}.$ (32)

This tool suggests three different types of dark energy behavior: quintessence
type $(\omega>-1)$ corresponding to its negative curvature (i.e. below the
$\Lambda CDM$ line), phantom type $(\omega<-1)$ relating to its positive
curvature (i.e. above the $\Lambda CDM$ line), and $Om(z)=\Lambda CDM$
corresponding to zero curvature. Fig. 10a reveal that the trajectories of
$Om(z)$ exhibit an increase as redshift $z$ decreases, indicating a negative
curvature pattern that confirms our model as a quintessence DE model for all
observational datasets. Our model collapses at 1 in the future i.e. at $z=-1$.

##  V  Thermodynamical behavior

This section examines the application of the generalized second law (GSL) of
thermodynamics to our model [89, 90]. It is important to note that
thermodynamics dictates that the entropy of isolated systems cannot decrease.
We investigate the computation of the total entropy $S$ for this model,
assuming that the apparent horizon is connected to temperature and entropy in
a way that is similar to the black hole event horizon. Consequently, the GSL
requires that the total entropy $S$ encompasses the entropy of all sources
[91]. As the Universe evolves, the rate of the entropy change for both the
fluid within it and the horizon must remain at zero or increase. In our
examination of recent and upcoming times, the total entropy is the combination
of the entropy contribution from matter $(S_{m})$ and the apparent horizon
$(S_{h})$, such that $S=S_{m}+S_{h}$, where $S_{h}=\frac{k_{B}A}{4l_{pl}^{2}}$
represents the entropy of the apparent horizon and $S_{m}$ represents the
entropy of pressure-less matter. The horizon’s area and Planck’s length are
denoted by $A$ and $l_{Pl}$, respectively, and the Boltzmann constant is
represented by $k_{B}$. The apparent horizon’s size can be determined using
the formula $A=4\pi r_{h}^{2}$, where
$r_{h}=\frac{1}{\sqrt{(}H^{2}+ka^{-2})}$. Since we’re focusing on a spatially
flat model $(k=0)$, we can assume that $r_{h}=H^{-1}$. This means that the
entropy of the horizon can be expressed as

$S_{h}=\frac{k_{B}\pi}{l_{pl}^{2}H^{2}}$ (33)

The Gibbs relation describes the entropy of matter within the dynamic apparent
horizon as [92, 93]

$TdS_{m}=d(\rho V)+pDV=Vd\rho+(p+\rho)dV$ (34)

The equation assumes that the temperature of the fluid, denoted by $T$, is
equal to the temperature at the horizon $(T_{h})$, which is calculated as
$T_{h}=\frac{1}{2\pi r_{h}}$, where $r_{h}$ is the radius of the horizon. $V$,
on the other hand, represents the spatial volume bounded by the horizon and is
equal to $\frac{4\pi}{3}r_{h}^{3}$. Differentiating Eqs. (33) and (34) and
adding them gives us the value of $\dot{S_{h}}+\dot{S_{m}}$.

(a)

(b)

Figure 10:  The plots depict the behavior of the Om diagnostic and the
variation of the temperature of the fluid in terms of the sum of the rates of
change of entropy for matter and the apparent horizon i.e.
T($\dot{S_{h}}+\dot{S_{m}}$).

$\dot{S_{h}}+\dot{S_{m}}$ remains non-negative throughout the entire history
of the universe. This implies that the sum of the rates of change of entropy
for matter and the horizon denoted as $\dot{S_{h}}$ and $\dot{S_{m}}$
respectively, is always greater than or equal to zero, as shown in Fig. 10b.
Therefore, the Generalized Second Law of thermodynamics is always satisfied
for a spatially flat FRW universe containing pressure-less dark matter within
the region enclosed by the apparent horizon.

##  VI Conclusion

We have examined a cosmological model of a homogeneous and isotropic universe
within the $f(R,T)$ gravity framework using a specific differential equation,
which is parameterized in terms of the energy density of the scalar field
$\rho_{\phi}$. By employing this parameterization, we have reconstructed the
Hubble parameter, $H(z)$, to analyze the evolution of the universe in a flat
FLRW space-time. The four model parameters are constrained using the Hubble
data, Pantheon data, BAO and their joint data. In order to comprehend cosmic
evolution, certain significant cosmological parameters have been formulated as
a function of redshift $z$ and illustrated through graphical representation.
Additionally, our obtained model has been compared with the standard model
$\Lambda$CDM.

The contour plots depict the best-fit values of the model parameters
$\Omega_{m0}$, $\Omega_{\phi 0}$, $H_{0},$ and $\alpha$ with $1-\sigma$ and
$2-\sigma$ errors for the Hubble, Pantheon, BAO, and their joint datasets (see
Figs. 1-4). The MCMC analysis indicates that the optimal values of the model
parameters derived from the Hubble dataset are slightly dissimilar to those
obtained from the Pantheon dataset. The best-fit values obtained from these
datasets align with the $\Lambda$CDM model. The error bar plots for both the
Hubble dataset and the SNeIa Union 2.1 compilation dataset show a good fit
when our model is correlated with $\Lambda$CDM in Fig. 5. Our model is an
ever-accelerating model and its present value $q_{0}\approx-1.004$ aligns with
the available observational data. The behavior of the jerk parameter is
consistent with the $\Lambda$CDM at present $z=0$ (see Fig. 6).

The evolution of the EoS parameter of the scalar field, $\omega_{\phi}$, seems
to be coinciding with the $\Lambda$CDM in the past and in the quintessence
region at present before collapsing in the future for a short period of time
and emerging to $\omega_{\phi}\to-1$ at $z=-1$ (See Fig. 7a). The behavior of
$\omega_{eff}$ is depicted in Fig. 7b, which can be observed that it remains
in the quintessence region in the past and at present before collapsing to $0$
at $z=-1$. Fig. 8 depicts the behavior of energy density, the energy density
of the scalar field, and pressure. $\rho_{eff}$ and $\rho_{\phi}$ decrease
with time and collapses to $0$ in future i.e at $z=-1$. The effective
pressure, $p_{eff}$, is negative at present ($z=0$) and attains a extremum
value before collapsing to $0$ in the late times ($z=-1$). Therefore, this
model represents the quintessence model at present.

The null, dominant, and strong energy conditions have been determined through
this analysis. In these conditions, the SEC is especially important as it
characterizes gravity’s attractive or repulsive behavior in the context of
modified gravity theory. Based on the observed values, it has been seen that
NEC and DEC conform to the conditions derived from the Raychaudhuri equations,
while SEC does not meet these conditions. The SEC violation indicates the
accelerated expansion of the universe (See Fig. 9). Om diagnostic tool is used
to compare our obtained model with the standard model $\Lambda$CDM. The
behavior of the model is depicted in Fig. 10a, with a comparison made between
the standard $\Lambda$CDM model and the SCDM model. This model shows a
quintessence model in late times because the variation of the Om diagnostic
curve having negative slope in the Fig. 10a.

Finally, we have examined the $GSL$ of Thermodynamics. The entropy shows a
positive value in the early stages of the evolution and a negative value in
the later stages before converging to $0$. Therefore, the $GSL$ holds good
during the initial stages of the evolution, but it becomes invalid during the
later stages of the evolution (See Fig. 10b). Thus, we conclude that our model
is an eternal accelerated expanding model according to the observational
analysis.

## References

  * [1] P. A. R. Ade et al. [BICEP2], Phys. Rev. Lett. 112, no.24, 241101 (2014).
  * [2] P. A. R. Ade et al. [Planck], Astron. Astrophys. 594, A13 (2016).
  * [3] N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020).
  * [4] A. G. Riess et al. [Supernova Search Team], Astron. J. 116, 1009-1038 (1998).
  * [5] A. G. Riess et al. [Supernova Search Team], Astrophys. J. 607, 665-687 (2004).
  * [6] A. G. Riess, L. G. Strolger, S. Casertano, H. C. Ferguson, B. Mobasher, B. Gold, P. J. Challis, A. V. Filippenko, S. Jha and W. Li, et al. Astrophys. J. 659, 98-121 (2007).
  * [7] S. Perlmutter et al. [Supernova Cosmology Project], Astrophys. J. 517, 565-586 (1999).
  * [8] D. H. Weinberg, M. J. Mortonson, D. J. Eisenstein, C. Hirata, A. G. Riess and E. Rozo, Phys. Rept. 530, 87-255 (2013).
  * [9] T. Padmanabhan, Phys. Rept. 380, 235-320 (2003).
  * [10] P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75, 559-606 (2003).
  * [11] E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753-1936 (2006).
  * [12] P. J. E. Peebles and B. Ratra, Astrophys. J. Lett. 325 (1988), L17.
  * [13] B. Ratra and P. J. E. Peebles, Phys. Rev. D 37 (1988), 3406.
  * [14] J. K. Singh, A. Singh, G. K. Goswami and J. Jena, Annals Phys. 443 (2022), 168958.
  * [15] J. K. Singh, Shaily, S. Ram, J. R. L. Santos and J. A. S. Fortunato, Int. J. Mod. Phys. D, 32, no. 7, 2350040 (2023).
  * [16] C. Armendariz-Picon, V. F. Mukhanov and P. J. Steinhardt, Phys. Rev. D 63, 103510 (2001).
  * [17] A. Sen, JHEP 07, 065 (2002).
  * [18] M. R. Setare, Eur. Phys. J. C 50, 991-998 (2007).
  * [19] N. Arkani-Hamed, H. C. Cheng, M. A. Luty and S. Mukohyama, JHEP 05, 074 (2004).
  * [20] F. Piazza and S. Tsujikawa, JCAP 07, 004 (2004).
  * [21] M. R. Setare, J. Sadeghi and A. R. Amani, Phys. Lett. B 660, 299-306 (2008).
  * [22] M. R. Setare, Phys. Lett. B 642, 421-425 (2006).
  * [23] A. Y. Kamenshchik, U. Moschella and V. Pasquier, Phys. Lett. B 511, 265-268 (2001).
  * [24] O. Bertolami, C. G. Boehmer, T. Harko and F. S. N. Lobo, Phys. Rev. D 75, 104016 (2007)
  * [25] O. Bertolami and J. Paramos, Phys. Rev. D 77, 084018 (2008)
  * [26] O. Bertolami, F. S. N. Lobo and J. Paramos, Phys. Rev. D 78, 064036 (2008)
  * [27] T. Harko, Phys. Lett. B 669, 376-379 (2008)
  * [28] T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82, 451-497 (2010).
  * [29] T. Clifton, P. G. Ferreira, A. Padilla and C. Skordis, Phys. Rept. 513, 1-189 (2012).
  * [30] A. De Felice and S. Tsujikawa, Living Rev. Rel. 13, 3 (2010).
  * [31] S. Nojiri and S. D. Odintsov, Phys. Rept. 505, 59-144 (2011).
  * [32] S. Nojiri and S. D. Odintsov, TSPU Bulletin N8(110), 7-19 (2011)
  * [33] S. A. Appleby, R. A. Battye and A. A. Starobinsky, JCAP 06, 005 (2010)
  * [34] P. Kanti, J. Rizos and K. Tamvakis, Phys. Rev. D 59, 083512 (1999)
  * [35] S. Nojiri, S.D. Odintsov, Phys. Rev. D 78, 046006 (2008)
  * [36] K. Bamba, S. Nojiri and S. D. Odintsov, JCAP 10, 045 (2008)
  * [37] S. Nojiri and S. D. Odintsov, Gen. Rel. Grav. 36, 1765-1780 (2004)
  * [38] S. Nojiri and S. D. Odintsov, eConf C0602061, 06 (2006)
  * [39] K. Nozari and T. Azizi, Phys. Lett. B 680, 205-211 (2009)
  * [40] B. C. Paul, P. S. Debnath and S. Ghose, Phys. Rev. D 79, 083534 (2009)
  * [41] T. Chiba, Phys. Lett. B 575, 1-3 (2003)
  * [42] A. D. Dolgov and M. Kawasaki, Phys. Lett. B 573, 1-4 (2003)
  * [43] L. Amendola, D. Polarski and S. Tsujikawa, Phys. Rev. Lett. 98, 131302 (2007)
  * [44] L. Amendola, R. Gannouji, D. Polarski and S. Tsujikawa, Phys. Rev. D 75, 083504 (2007)
  * [45] A. W. Brookfield, C. van de Bruck and L. M. H. Hall, Phys. Rev. D 74, 064028 (2006)
  * [46] T. Harko, F. S. N. Lobo, S. Nojiri and S. D. Odintsov, Phys. Rev. D 84, 024020 (2011).
  * [47] M. J. S. Houndjo, F. G. Alvarenga, M. E. Rodrigues and D. F. Jardim, Eur. Phys. J. Plus 129, 171 (2014)
  * [48] M. Sharif and M. Zubair, JCAP 03, 028 (2012) [erratum: JCAP 05, E01 (2012)]
  * [49] M. Sharif and M. Zubair, J. Exp. Theor. Phys. 117, 248-257 (2013)
  * [50] F. G. Alvarenga, M. J. S. Houndjo, A. V. Monwanou and J. B. C. Orou, J. Mod. Phys. 4, 130-139 (2013)
  * [51] M. Sharif and M. Zubair, JHEP 12, 079 (2013).
  * [52] H. Shabani and M. Farhoudi, Phys. Rev. D 88, 044048 (2013)
  * [53] M. Kiran and D. R. K. Reddy, Astrophys. Space Sci. 346, 521-524 (2013)
  * [54] M. Sharif and M. Zubair, Astrophys. Space Sci. 349, 457-465 (2014)
  * [55] M. J. S. Houndjo, Int. J. Mod. Phys. D 21, 1250003 (2012).
  * [56] D. R. K. Reddy, R. L. Naidu, K. Dasu Naidu and T. Ram Prasad, Astrophys. Space Sci. 346, 261-265 (2013)
  * [57] E. H. Baffou, A. V. Kpadonou, M. E. Rodrigues, M. J. S. Houndjo and J. Tossa, Astrophys. Space Sci. 356, no.1, 173-180 (2015)
  * [58] P. H. R. S. Moraes, Astrophys. Space Sci. 352, 273-279 (2014)
  * [59] S. Nojiri, S. D. Odintsov and V. K. Oikonomou, Phys. Rept. 692, 1-104 (2017).
  * [60] M. J. S. Houndjo and O. F. Piattella, Int. J. Mod. Phys. D 21, 1250024 (2012).
  * [61] M. Jamil, D. Momeni, M. Raza and R. Myrzakulov, Eur. Phys. J. C 72 (2012), 1999.
  * [62] F. G. Alvarenga, A. de la Cruz-Dombriz, M. J. S. Houndjo, M. E. Rodrigues and D. Sáez-Gómez, Phys. Rev. D 87, no.10, 103526 (2013) [erratum: Phys. Rev. D 87, no.12, 129905 (2013)].
  * [63] Z. Haghani, T. Harko, F. S. N. Lobo, H. R. Sepangi and S. Shahidi, Phys. Rev. D 88, no.4, 044023 (2013).
  * [64] S. D. Odintsov and D. Sáez-Gómez, Phys. Lett. B 725, 437-444 (2013).
  * [65] H. Shabani and M. Farhoudi, Phys. Rev. D 90, no.4, 044031 (2014).
  * [66] J. K. Singh, K. Bamba, R. Nagpal and S. K. J. Pacif, Phys. Rev. D 97, no.12, 123536 (2018).
  * [67] H. Shabani, Int. J. Mod. Phys. D 26 (2017) no.10, 1750120.
  * [68] J. K. Singh, H. Balhara, K. Bamba and J. Jena, JHEP 03, 191 (2023) [erratum: JHEP 04 (2023), 049].
  * [69] L. D. Landau and E. M. Lifschits, Pergamon Press, 1975, ISBN 978-0-08-018176-9
  * [70] J. K. Singh and R. Nagpal, Eur. Phys. J. C 80, no.4, 295 (2020).
  * [71] D. M. Scolnic et al. [Pan-STARRS1], Astrophys. J. 859, no.2, 101 (2018).
  * [72] Shaily, M. Zeyauddin and J. K. Singh, [arXiv:2207.05076 [gr-qc]].
  * [73] C. Blake, E. Kazin, F. Beutler, T. Davis, D. Parkinson, S. Brough, M. Colless, C. Contreras, W. Couch and S. Croom, et al. Mon. Not. Roy. Astron. Soc. 418 (2011), 1707-1724
  * [74] W. J. Percival et al. [SDSS], Mon. Not. Roy. Astron. Soc. 401 (2010), 2148-2168
  * [75] D. J. Eisenstein et al. [SDSS], Astrophys. J. 633 (2005), 560-574
  * [76] R. Giostri et al., J. Cosmol. Astropart. Phys. 2012(03) (2012) 027
  * [77] Y. L. Bolotin, V. A. Cherkaskiy, O. A. Lemets, D. A. Yerokhin and L. G. Zazunov, [arXiv:1502.00811 [gr-qc]].
  * [78] R. D. Blandford, M. A. Amin, E. A. Baltz, K. Mandel and P. J. Marshall, ASP Conf. Ser. 339, 27 (2005).
  * [79] D. Rapetti, S. W. Allen, M. A. Amin and R. D. Blandford, Mon. Not. Roy. Astron. Soc. 375, 1510-1520 (2007).
  * [80] T. Chiba and T. Nakamura, Prog. Theor. Phys. 100, 1077-1082 (1998).
  * [81] M. Visser, Class. Quant. Grav. 21, 2603-2616 (2004).
  * [82] O. Luongo, Mod. Phys. Lett. A 28, 1350080 (2013).
  * [83] A. Al Mamon and K. Bamba, Eur. Phys. J. C 78, no.10, 862 (2018).
  * [84] U. Alam, V. Sahni, T. D. Saini and A. A. Starobinsky, Mon. Not. Roy. Astron. Soc. 344, 1057 (2003).
  * [85] V. Sahni, T. D. Saini, A. A. Starobinsky and U. Alam, JETP Lett. 77, 201-206 (2003).
  * [86] Y. G. Gong and A. Wang, Phys. Rev. D 75, 043520 (2007).
  * [87] V. Sahni, A. Shafieloo and A. A. Starobinsky, Phys. Rev. D 78, 103502 (2008).
  * [88] C. Zunckel and C. Clarkson, Phys. Rev. Lett. 101, 181301 (2008).
  * [89] A. Sheykhi and B. Wang, Phys. Lett. B 678, 434-437 (2009).
  * [90] M. R. Setare and E. C. Vagenas, Phys. Lett. B 666, 111-115 (2008).
  * [91] J. D. Barrow, Nucl. Phys. B 310, 743-763 (1988).
  * [92] E. K. Li, Y. Zhang, J. L. Geng and P. F. Duan, Gen. Rel. Grav. 47, no.11, 136 (2015).
  * [93] G. Izquierdo and D. Pavon, Phys. Lett. B 633, 420-426 (2006).


# Incorporating baseline covariates to validate surrogate endpoints with a
constant biomarker under control arm

Emily Roberts1∗, Michael Elliott1,2, Jeremy M. G. Taylor1
1Department of Biostatistics, University Michigan
1415 Washington Heights Ann Arbor, MI 48109
<EMAIL_ADDRESS>
2 Survey Methodology Program, Institute for Social Research Ann Arbor, MI
48109

### Abstract

A surrogate endpoint $S$ in a clinical trial is an outcome that may be
measured earlier or more easily than the true outcome of interest $T$. In this
work, we extend causal inference approaches to validate such a surrogate using
potential outcomes. The causal association paradigm assesses the relationship
of the treatment effect on the surrogate with the treatment effect on the true
endpoint. Using the principal surrogacy criteria, we utilize the joint
conditional distribution of the potential outcomes $T$, given the potential
outcomes $S$. In particular, our setting of interest allows us to assume the
surrogate under the placebo, $S(0)$, is zero-valued, and we incorporate
baseline covariates in the setting of normally-distributed endpoints. We
develop Bayesian methods to incorporate conditional independence and other
modeling assumptions and explore their impact on the assessment of surrogacy.
We demonstrate our approach via simulation and data that mimics an ongoing
study of a muscular dystrophy gene therapy.

### Keywords

Bayesian methods, principal stratification, subgroup effects, surrogate
endpoints

##  1 Introduction

Although randomized clinical trials are largely considered the gold standard
to evaluate treatment efficacy, methods that lower trial cost and shorten the
length of the study are often sought after in the medical field. Surrogate
endpoints $S$ are biologically plausible intermediate outcomes that are
strongly related to the true outcome of interest $T$. In a trial, these
endpoints may be measured earlier or more effectively to quickly disperse
treatments to patients. Popular examples of potential surrogate endpoints
include CD4 blood counts for HIV mortality and immune responses for vaccine
efficacy. It is crucial to collect data and validate such an endpoint before
using in a large-scale trial.

Prentice’s landmark 1989 paper proposed criteria in a single trial setting to
evaluate statistical surrogates: that the surrogate $S$ both be related to the
outcome and that it captures the effect of the treatment $Z$ on $T$.1 Other
criteria have since been proposed, as it has been shown simple criteria may
not ensure a seemingly useful surrogate will predict a beneficial treatment
effect.2 $S(z)$ and $T(z)$ refer to the endpoint values had the treatment,
possibly counter-factually, been assigned to level $z$. Since $S$ is measured
after treatment assignment, conditioning on the surrogate distorts the causal
pathway and interpretation of the treatment effect in a regression model.
Using Rubin’s potential outcome causal framework,3,4 principal surrogacy
proposed the solution of using both potential intermediate outcomes by
considering the surrogate values under each treatment as pre-treatment
variables.5 Since both surrogate outcomes are hypothetically determined prior
to randomization, surrogates can be evaluated based on principal causal
effects defined on the distribution of $T(1)-T(0)$ conditional on principal
strata defined by $S(1)-S(0)$. Gilbert and Hudgens (forward as GH6) elaborated
on causal necessity, as proposed by Frangakis and Rubin, by terming average
causal necessity and average causal sufficiency. These require that there be
no average effect of the treatment on the true outcome in the strata where
there is no average effect on the surrogate, and similarly that there exists
an average treatment effect on the true outcome in the strata where is there
an average effect on the surrogate. Further, they define the Causal Effect
Predictiveness (CEP) curve as a visualization of both requirements across the
value of the surrogate endpoints $(s_{1},s_{0})$.

Our work is motivated by an ongoing study of a muscular dystrophy treatment.7
In this trial of a new gene transfer therapy, patients received a micro-
dystrophin transgene to produce the micro-dystrophin protein. The potential
surrogate $S$ is micro-dystrophin expression as measured by western blot
methods, and the primary outcome of interest $T$ is the North Star Ambulatory
Assessment (NSAA) functional score on a continuous scale. $S$ is measured at
only one time point, while $T$ is measured before randomization as well as
after the gene transfer therapy. Since muscular growth and deterioration due
to the disease have major impact on physical movement during childhood, both
baseline ambulatory ability and age are important to take into consideration.
In this setting, patients do not produce significant amounts of micro-
dystrophin protein at baseline, so it can be assumed that the value of the
surrogate under placebo is approximately equal to 0. This scenario where
$S(0)$ is fixed to 0 also commonly arises in vaccine efficacy studies. Since
those in the placebo group necessarily have no immune response without the
vaccine antigens, GH refer to this simplified setting as the constant
biomarker case.

Quantities related to vaccine efficacy were developed in the HIV and pertussis
settings,8,9,10 and GH were among the first to formalize the surrogate
validation methodology in a HIV vaccine efficacy trial. Still, a major
challenge of characterizing these causal effect summaries is dealing with
nonidentified parameters arising from use of potential outcomes, so Follmann12
suggested the closeout placebo vaccination design to avoid the unobserved
outcomes. GH focused on modeling assumptions to identify the causal quantities
and generalized previous work of the baseline immunogenicity predictor (BIP)
$W$ to estimate the missing $S(z)$ value. Related work has proposed augmented
trial design ideas such as the baseline surrogate measure and the cross-over
design, respectively, and other authors have imposed conditional independence
assumptions of the outcomes13. Subsequently, several authors have addressed
particular models and designs for vaccines and immune correlates of
protection.14,15,16,17,18

Under the general Bayesian paradigm, with specific assumptions about parameter
values, both Zigler and Belin19 and Conlon, Taylor, and Elliott20 proposed to
consider the full joint distribution of potential outcomes to create the CEP
curve using imputation strategies applicable to settings beyond the constant
biomarker case. Work in the frequentist setting by Alonso, Van der Elst, and
Meyvisch utilized potential outcomes and the information-theoretic framework
to propose a surrogate predictive function with a two-step procedure for
dealing with non-identifiability.21 Making the constant biomarker assumption
results in fewer missing potential outcomes to impute, which allows us to
focus on the sensitivity of modeling assumptions. In this work, we build upon
previous models of the joint distribution of potential outcomes $S(0),\ S(1),\
T(0),\ T(1)$ by exploring trial and modeling considerations when controlling
only three of these endpoints as applicable to our motivating clinical trial.
20,22

We propose techniques and design considerations with the goal of achieving
gains in estimation efficiency. Since the surrogate and true outcome values
are only observed for the assigned treatment, we consider the counterfactual
outcomes as missing data and implement an imputation strategy for estimation.
We compare this algorithm to instead using only observed data and prior
distributions for nonidentified correlation parameters. The novelty of this
work is the incorporation of baseline covariates with two objectives in mind:
first, conditioning on baseline covariates may improve the plausibility of
conditional independence assumptions, and second, it allows us to make
inference about whether there are subgroups of the population for whom the
quality of the surrogate varies. For the latter, we propose to stratify the
previously marginal estimands for validation by conditioning on patient
characteristics. In our application, we focus on one such example of a
baseline covariate, namely that of the true outcome measured pre-treatment,
which is similar to the BSM proposed by other authors. This particular
measurement allows for multiple definitions of the true outcome of interest
and other trial design decisions. We recognize that by viewing the baseline as
a noisy estimate of $T(0)$ similar to a measurement error problem, it could
provide improved identifiability or yield more informative prior distributions
for nonidentified parameters.

In Section 2, we propose the model and incorporation of baseline covariates in
the surrogate setting. We define the conditional surrogacy validation metrics
and suggest potential reasons to use the covariates such as to make
conditional independence assumptions and raise consideration for how to define
the trial endpoint. In Section 3, we describe the proposed Bayesian estimation
methods using either an imputation scheme or observed data algorithm.
Simulation studies are shown in Section 4, and our data example is explored in
Section 5 before a concluding discussion in Section 6.

##  2 The Model

Using the causal association framework, we first consider the joint
distribution of three continuous potential outcomes under a binary treatment
$Z$. Since $S(0)=0$, we assume a multivariate normal distribution of the
counterfactual surrogate and true outcomes for each subject:

$\left(\begin{array}[]{c}S(1)\\\ T(0)\\\ T(1)\end{array}\right)\sim
MVN\left(\left(\begin{array}[]{c}\delta_{1}\\\ \delta_{2}\\\
\delta_{3}\end{array}\right),\left(\begin{array}[]{ccc}\sigma_{S1}^{2}&\rho_{10}\sigma_{S1}\sigma_{T0}&\rho_{11}\sigma_{S1}\sigma_{T1}\\\
&\sigma_{T0}^{2}&\rho_{T}\sigma_{T0}\sigma_{T1}\\\ &&\sigma_{T1}^{2}\\\
\end{array}\right)\right)$ (1)

###  2.1 Assumptions

In this setting, we focus on clinical trial scenarios where we can assume
$S(0)=0$ (the methods can be extended to more general settings). The causal
inference assumptions we make are the Stable Unit Treatment Values Assumption
and ignorable treatment assignment, meaning the potential outcomes for any
individual do not vary with the treatments assigned to other individuals, and
treatment assignment is independent of potential outcomes conditional on all
covariates ($P(Z=1|T(0),T(1),X)=P(Z=1|X)$), respectively. Since we do not
observe combinations of joint outcomes $\\{T(0),S(1)\\}$ or $\\{T(0),T(1)\\},$
the correlation parameters $\rho_{10},\rho_{T}$ are not identified. We will
consider various approaches to obtain identifiability through the use of
proper priors, conditional independence assumptions, and/or fixing unknown
parameters via sensitivity analyses, and consider how our models may be
adapted if additional baseline data is available. Then we use the specified
joint model to impute the missing counterfactual values. This potential
outcomes approach captures the causal associations for validation.

The multivariate normality assumption provides many convenient results by
being analytically tractable and allowing for closed form quantities for
surrogate validation which we describe below. As this may not hold in
practice, we later verify the sensitivity of this distributional assumption
and assess the robustness of the results in the presence of model
misspecification. Other work has incorporated copula models for non-normal
data.23

###  2.2 Baseline Covariates

Some estimates for surrogacy quality have wide confidence bands that make
definitive recommendations difficult. While most causal metrics are reported
marginally, it may be beneficial to use baseline covariates $X$ in the
analysis. One interest is to assess effect modification: that is, if there
exist subgroups of patients to determine for whom the surrogate will work
particularly well for the true outcome, such as for males or those who are
young. This has the potential to reduce the risk of observing the surrogate
fail in a certain patient population after approval for use in subsequent
trials. It is also possible that covariates would help predict membership of
principal strata. $X$ may explain dependence or confounding between $S$ and
$T$ that may occur in finite samples, even after trial randomization.
Statistical benefits may be seen in the estimation accuracy as well via more
accurate imputation of the missing counterfactual values to both reduce bias
and gain efficiency. Finally, conditioning on $X$ might be expected to reduce
correlations amongst the potential outcomes, allowing us to make stronger
conditional independence statements after conditioning on baseline covariates.
Note that these latter examples’ effects can be reported as either conditional
on $X$, or integrated over the empirical distributions of covariates to
provide marginal estimates.

The conditional model can be written with effects of $X$ in the mean structure
(therefore the parameters in this covariance structure differ and $\theta$
represents the conditional correlations)

$\left(\begin{array}[]{c|c}S(1)&\\\ T(0)&X\\\ T(1)\end{array}\right)\sim
N\left(\left(\begin{array}[]{c}\omega_{1}+\omega_{2}X\\\
\omega_{3}+\omega_{4}X\\\
\omega_{5}+\omega_{6}X\end{array}\right),\left(\begin{array}[]{ccc}\epsilon_{S1}^{2}&\theta_{10}\epsilon_{S1}\epsilon_{T0}&\theta_{11}\epsilon_{S1}\epsilon_{T1}\\\
&\epsilon_{T0}^{2}&\theta_{T}\epsilon_{T0}\epsilon_{T1}\\\
&&\epsilon_{3}^{2}\\\ \end{array}\right)\right)$ (2)

We note that this model still has two nonidentified parameters, $\theta_{10}$
and $\theta_{T}$. In the above model $X$ is a scalar, but it could be
generalized to a vector. Furthermore, one of the $X$ components might be known
to be highly related to either $S$ or $T$ (such as a pre-treatment
measurement). In these specific cases that provide additional information, it
may be feasible to make further assumptions about the model structure.

###  2.3 Surrogacy Validation

The validation causal quantities derived from conditioning on strata of the
surrogate can be written as a function of the model parameters. These
quantities can be viewed graphically in the causal effect predictiveness (CEP)
surface as a line with intercept and slope based on causal effects as the
difference in surrogate potential outcomes $S(1)-S(0)=s$ on the $x$-axis and
difference in the expected, conditional true outcomes
$E(T(1)-T(0)|S(1)-S(0)=s)$ on the $y$-axis. In the case of Gaussian
distributions, as in equation 1, $E(T(1)-T(0)|S(1)-S(0)=s)$ is linear in $s$
and has the form $=\gamma_{0}+\gamma_{1}s$. By displaying expected change in
potential outcomes, conditional on the actual surrogate change $s$, the plots
demonstrate if the surrogate is valid, meaning small (large) causal effects on
a surrogate are associated with small (large) causal effects on the outcome.
When the distribution of outcomes is multivariate normal, average causal
necessity and average causal sufficiency are fulfilled if $\gamma_{0}=0$, the
expected change in true outcome when there is no change in the surrogate
outcome at the origin, and $\gamma_{1}\neq 0$, the expected change in true
outcome when there is a nonzero change in the surrogate outcome. Under the
multivariate normal distribution in equation 1, these values from the
conditional expectation can be written as

$\gamma_{0}=(\mu_{T1}-\mu_{T0})-\gamma_{1}\mu_{S1}=(\delta_{3}-\delta_{2})-\gamma_{1}\delta_{1}\
\ \ \ \ \ \ \ \
\gamma_{1}=\frac{\rho_{11}\sigma_{T1}-\rho_{10}\sigma_{T0}}{\sigma_{S1}}$

or when incorporating baseline covariates,
$\gamma_{0}=(\mu_{T1|X}-\mu_{T0|X})-\gamma_{1}\mu_{S1|X}$ where
$\gamma_{1}=\frac{\theta_{11}\epsilon_{T1}-\theta_{10}\epsilon_{T0}}{\epsilon_{S1}}$.
Our goal is to estimate all parameters in the distribution so we can calculate
$\gamma_{0}$ and $\gamma_{1}$ and determine if $S$ is a valid surrogate for
$T$. We can understand $\gamma_{1}$ by rewriting the quantity as
$\frac{\omega_{11}-\omega_{10}}{\sigma_{S1}}$ and can consider when
$\omega_{11}>\omega_{10}$ that the slope will be positive, and when the ratio
$\omega_{11}/\omega_{10}$ is larger, the magnitude of $\gamma_{1}$ is larger.
Another way to look at this term is as
$\frac{Cov(T(1),S(1))-Cov(T(0),S(1))}{Var(S(1))}$, where the sign and
magnitude of $\gamma_{1}$ is determined by the covariances between $S(1)$ and
the true outcomes.

For scenarios where we believe the surrogate works particularly poorly for
certain patient groups, we would be interested in a stratified CEP curve. In
other settings, we may simply incorporate covariates as an intermediate step
while remaining interested in the marginal CEP curve. To denote the
difference, let $\gamma_{0,C}$ and $\gamma_{1,C}$ correspond to a model fit
using $X$, compared to the marginalized estimates that is accomplished by
empirically averaging over the distribution of $X$. Later, we will further
differentiate these respective conditional $C$ and marginal $M$ models based
on how we define the outcome. From equation 1, $T(1)-T(0)|S(1),x$ has a normal
distribution, so the expected value will be written as
$\gamma_{0,C}+\gamma_{1,C}s$. In this expression $\gamma_{0,C}$ can depend on
the covariate $X$, but by the assumptions of the model $\gamma_{1,C}$ (i.e.
the covariance) does not depend on $X$. From these values, the marginalization
is written

$\int_{x}E(T(1)-T(0)|S(1),X=x)f(X|S(1)=s)dx=\int_{x}\frac{E(T(1)-T(0)|S(1),x)f(S(1)|x)f(X)}{f(S(1))}dx$
(3)

using Bayes rule (see Web Appendix A). Once we obtain these quantities, we
plot the marginal effect over values of $s$ and summarize it by fitting a
linear model to estimate $\gamma_{0,M}$ and $\gamma_{1,M}$.

###  2.4 Subgroups in Surrogacy and Treatment Effects

Based on our definition of $\gamma_{0},$ our concept of surrogacy subgroups is
analogous to heterogeneous treatment effects existing. Since the quantity
$\gamma_{0,C}$ depends on $x$, a surrogacy subgroup effect will occur only
when there is an interaction with the treatment effect on either the surrogate
or the true outcome. In contrast, we are not assuming the absence of
unobserved heterogeneity in the sense of a sharp null existing.24

###  2.5 Conditional Independence

Conditional independence assumptions are frequently made in causal inference
and for surrogate endpoint validation in particular (see Conlon et al., for
examples of assumptions made)18,25,26,6. Briefly, when $S(0)$ is not
restricted to the value of 0, common assumptions among the endpoints include
strong statements about the outcomes such as 1) $T(0)\perp T(1)|S(0),S(1)$; 2)
$S(0)\perp T(1)|S(1)$; 3) $S(1)\perp T(0)|S(0)$; and intuitively weaker 4)
$S(0)\perp T(1)|S(1),T(0)$ and 5) $S(1)\perp T(0)|S(0),T(1)$.

In our setting, the reasonable assumption that is still feasible is a
collapsed version of 5: $S(1)\perp T(0)|T(1)$. Determining if a conditional
independence assumption is plausible is context-dependent. Examining the
assumption, it is not unreasonable to believe that given the true outcome
under $Z=1$, the surrogate for $Z=1$ is independent of the true outcome under
the opposing treatment. This assumption implies the relationship among
equation 2 parameters that $\theta_{T}=\theta_{10}/\theta_{11}$. In the
multivariate normal setting, this combination of correlation parameters on the
joint normal scale is derived by setting a conditional covariance term to $0$
(described in Web Appendix A), which corresponds to a zero element in the
inverse of the correlation matrix. After adjusting for baseline covariates
$X$, conditional independence of the outcomes may be more likely to hold,
which in turn reduces the number of parameters to estimate by one (a
restatement of equation 2 in which $\theta_{10}$ has been replaced by
$\theta_{11}\times\theta_{T}$):

$\left(\begin{array}[]{c|c}S(1)&\\\ T(0)&X\\\ T(1)\end{array}\right)\sim
N\left(\left(\begin{array}[]{c}\omega_{1}+\omega_{2}X\\\
\omega_{3}+\omega_{4}X\\\
\omega_{5}+\omega_{6}X\end{array}\right),\left(\begin{array}[]{ccc}\epsilon_{S1}^{2}&\theta_{11}\theta_{T}\epsilon_{S1}\epsilon_{T0}&\theta_{11}\epsilon_{S1}\epsilon_{T1}\\\
&\epsilon_{T0}^{2}&\theta_{T}\epsilon_{T0}\epsilon_{T1}\\\
&&\epsilon_{T1}^{2}\\\ \end{array}\right)\right)$ (4)

This potentially increases efficiency and helps with identifiability since
there is only one nonidentified parameter ($\theta_{T}$).

###  2.6 Design Considerations and Defining the Endpoint

In our motivating clinical trial, one of our baseline covariates is actually a
pre-treatment measurement of the outcome $T.$ Because of this, we could choose
to define the true outcomes several ways with the possible benefit of
maximizing efficiency. For example, we could define the outcomes $T(0),T(1)$
as the original values of the endpoints and use $X$ for subgroup analysis.
Alternatively, we could analyze the outcomes as the change from baseline
measurement $X$ to later measurements, denoted as $T^{D}(0)$ and $T^{D}(1)$.
We would like to know if it is advantageous to define a trial outcome in terms
of change from baseline compared to using the baseline value as a covariate.
To be explicit about the quantities we are estimating, we can outline the
relevant joint and conditional distributions based on the choice of endpoint
definition and involvement of $X$. All four methods we consider are based on
the joint 4x4 distribution from equation 1 extended to include $X$. In the
special setting where our covariate is a pre-measurement of $T(0)$, we briefly
consider other strong assumptions that could be made. For example, in the
strongest case we could plug in $X$ in for $T(0)$ and gain identifiability.
Alternatively, we could view $X$ and $T(0)$ as repeated measures and assume
certain parameters are equal, such as identifiable means and variances or
nonidentified correlations. We consider how these relate to conditional
independence assumptions proposed by other authors in future work.

In the setting we are considering, there is a pre-planned analysis of the
final clinical trial data, and defining the endpoint is necessary to test for
a treatment effect. Like we have suggested for assessing surrogacy, this step
involves fitting marginal or conditional models with the original outcome or
difference from baseline endpoint. For efficiency or to incorporate treatment
effect subgroups, we may choose to condition on baseline covariates $X$. We
enumerate the potential analysis models for the treatment effect based on the
observed outcome $T$ and the corresponding surrogate validation models for the
potential outcomes $T(0)$ and $T(1)$ conditional on $S(1)=s$. These are
differentiated by the endpoint:

  1. 1.

Original outcome: $T(0),T(1)$ Treatment effect model
$T_{i}=\beta_{0}+\beta_{1}Z_{i}+\epsilon_{i}$

Surrogate validation metric $E(T(1)-T(0)|s)=\gamma_{0,M}+\gamma_{1,M}s$

  2. 2.

Also condition on $X$: $T(0),T(1)|X$ Treatment effect model
$T_{i}=\beta_{2}+\beta_{3}Z_{i}+\beta_{4}X_{i}+\epsilon_{i}$

Surrogate validation metric $E(T(1)-T(0)|X,s)=\gamma_{0,C}+\gamma_{1,C}s$

  3. 3.

Difference from baseline: $T^{D}(0),T^{D}(1)$ Treatment effect model
$T^{D}_{i}=\beta_{5}+\beta_{6}Z_{i}+\epsilon_{i}$

Surrogate validation metric
$E(T^{D}(1)-T^{D}(0)|s)=\gamma_{0,M}+\gamma_{1,M}s$

  4. 4.

Also condition on $X$: $T^{D}(0),T^{D}(1)|X$ Treatment effect model
$T^{D}_{i}=\beta_{7}+\beta_{8}Z_{i}+\beta_{9}X_{i}+\epsilon_{i}$

Surrogate validation metric
$E(T^{D}(1)-T^{D}(0)|X,s)=\gamma_{0,C}+\gamma_{1,C}s$

Parameters $\beta_{1},\beta_{3},\beta_{6},$ and $\beta_{8}$ estimate the
respective treatment effects. We expect the first and third methods to produce
estimates of the same population, marginal treatment effect regardless of
subtracting off the baseline measurement, though finite sample equality is
unlikely to hold for even a randomized trial. These same considerations extend
to the surrogate validation framework, where some of these methods will
estimate the same marginal validation estimates. We will later consider a
further reason to thoughtfully define the endpoint, which is to determine
which scale it is reasonable to assume conditional independence.

##  3 Bayesian Methods

###  3.1 Imputation-Estimation Algorithms

While we observe only $S(0),T(0)$ for $n_{0}$ subjects and $S(1),T(1)$ for
$n_{1}$ subjects, our validation quantities involve correlations that can only
be calculated from counterfactual outcomes. In order to simultaneously
estimate the model parameters, address nonidentified terms, and to
appropriately propagate the uncertainty of imputing missing outcomes, we use a
Bayesian method for estimation. There are three types of variables that will
be iteratively drawn in the MCMC algorithm: the correlation parameters in $R$,
the model mean and variance parameters ($\mu$’s and $\sigma$’s), and the
missing potential outcomes. The imputation strategy is a full process that
iteratively imputes the missing potential outcomes and uses the posterior
distribution to draw values of the parameters.

We assume vague, normal priors for the identified mean parameters. Rather than
use an Inverse-Wishart prior or another method that would sample the entire
matrix at once that lacks needed flexibility, we implement a separation method
on the covariance matrix $\Sigma$. This decomposes the matrix $\Sigma=QRQ$
into standard deviation and correlation matrices where $R$ is a correlation
matrix with $1$’s on the diagonal to easily place less informative priors on
the identified terms.27 To ensure iterative draws satisfy the positive
definite constraint on $\Sigma$, the method uses the griddy Gibbs algorithm to
draw from the appropriate bounded posterior.20 Specifically, we compute the
posterior of each parameter over a set of realizable grid points and re-
evaluate over a region of high posterior density with more precision (finer
grid points) before randomly drawing the value for that iteration. We consider
different priors for the correlation parameter and find that it is important
to carefully choose on which parameters to place priors when implementing the
conditional independence constraint. Since we are in a setting with
nonidentified correlation parameters where the data will provide no direct
information about the true values of these parameters, we are careful to not
impose unreasonable prior distributions as we expect the posterior to mimic
the prior. We consider both vague and more informative Uniform and Beta priors
on the correlation terms, though the marginal distribution of each correlation
under positive definite constraints can become less straightforward. We assume
$S(1)\perp T(0)|T(1)$ (or its equivalent based on the exact model fit and
incorporation of $X$), by drawing suitable values of $\theta_{T}$ and
$\theta_{11}$ using the grid search. Essentially, as demonstrated in equation
4, the term $\theta_{10}$ is no longer involved in the likelihood, and the
product $\theta_{T}\times\theta_{11}$ takes its place.

###  3.2 Observed Data Algorithm

An alternative way to estimate the identifiable parameters is by using the
observed data likelihood and devising an MCMC algorithm to obtain draws from
the posterior distribution. This approach avoids imputing the counterfactual
values of $S$ and $T$. Since the data contain no information about
$\theta_{T}$ or $\theta_{10}$, we expect the posterior to match the
distribution of the chosen prior, provided the priors are independent. Thus,
we propose to draw the nonidentified correlation parameters directly from the
prior solely for the purposes of estimating $\gamma_{0}$ and $\gamma_{1}$.
Since the value of $\theta_{11}$ is estimable from the observed data, we use
the posterior distribution for its draw, but it is the only correlation
parameter drawn from a conditional distribution. When assuming conditional
independence, we replace the term $\theta_{10}$ with
$\theta_{T}\times\theta_{11}$ (after $\theta_{11}$ is drawn from the posterior
and $\theta_{T}$ is drawn directly from
Uniform(-1, 1)) in the same way as explained above when calculating values of
$\gamma_{1}$ and ensuring the matrix is positive definite. We then fit the
regression models on the observed data only to estimate the other parameters.

Here we carry out basic Bayesian estimation of the identified parameters for
comparability of uncertainty estimates to isolate the effect of using the
imputation scheme and priors for nonidentified parameters. We find that we can
bypass the full MCMC scheme intended to provide parameter estimates while
addressing the nonidentified parameters, and instead we can draw these
independently from the prior. This is related to work by Gustafson that
demonstrates any difference between the prior and posterior distribution for
non-identified parameters is due to prior dependence in the parameters.28
Using his transparent reparameterization here, there is no indirect learning
of the correlations outside of the positive definite and conditional
independence assumption constraints, which are still enforced with this
algorithm. In this setting, we expect results from this method to be generally
equivalent to imputation while being less computationally expensive. We also
note that using an MCMC scheme to estimate the variance and mean parameters
may not be necessary at all, and a maximization of the posterior or maximum
likelihood method may be used instead.

##  4 Simulation Studies

We explore the impact of using a baseline covariate in terms of efficiency and
making conditional independence assumptions. In particular, our simulations
are meant to assess how we define the true endpoint based on different
relationships with the baseline covariate. Using simulation studies, we
generate data that mimics a randomized trial, meaning for half of the subjects
assigned $Z=0$, we observe only $T(0),X$, and for the other half we observe
$S(1),T(1),X$. The five sets of generative parameters for equation 2 are shown
in Table 1. The five settings have different generating parameters that vary
the treatment effect and quality of the surrogate endpoint to demonstrate the
method’s performance over a variety of scenarios. Based on the model in
equation 2 for observed data of sample size $n=100$, we generate data that
fulfills conditional independence in all scenarios using the original endpoint
and conditional on $X$ setting. [Table 1]

The surrogate is valid marginally only for settings $A,B$. The covariate $X$
is normally distributed in settings $A$-$D$ and binary in $E$; settings $D,E$
represent the existence of a subgroup effect. After generating the data, we
fit the models described below and vary which conditional independence
assumption is made, if any. To perform surrogacy validation, we fit four
models of tri-variate normal distributions derived from the distribution of
$S(1),T(0),T(1),X$ (see Web Appendix B for details). The four analysis designs
are based on the distribution of the three outcomes and baseline covariates,
and the different parameterizations can be easily equated algebraically. For
any $X$ (we simulate $X\sim N(\delta_{4},\sigma_{X}^{2})$ or $X\sim
Bernoulli(0.5)$),

1. $\left(\begin{array}[]{c}S(1)\\\ T(0)\\\ T(1)\end{array}\right)$ 2. $\left(\begin{array}[]{c|c}S(1)&\\\ T(0)&X\\\ T(1)&\end{array}\right)$ 3. $\left(\begin{array}[]{c}S(1)\\\ T^{D}(0)\\\ T^{D}(1)\end{array}\right)$ 4. $\left(\begin{array}[]{c|c}S(1)&\\\ T^{D}(0)&X\\\ T^{D}(1)&\end{array}\right)$

We consider that conditional on baseline covariates based on context, we may
decide to introduce strong modeling assumptions. For example, when $X$ is a
pre-treatment measurement of $T$ (like a measurement-error prone value of
$T(0)$), priors may be informed by estimates of the observed correlation of
$X$ and $T(1)$, or we may be able to estimate some nonidentified correlations.
Currently we use a Beta prior (truncated between -0.4 to 1 with a positive
mean equal to 0.23) on the correlation between either $T(0),T(1)$ or
$T^{D}(0),T^{D}(1)$ and a Uniform(-1, 1) prior on the correlation between
$S(1),T(1)$ or $S(1),T^{D}(1)$ when conditional independence is not assumed.
To perform sensitivity analyses, we will both vary the prior distributions on
the nonidentified parameters to integrate over the range of plausible values
and fix the correlations to see at what boundaries the conclusions change.

These models allow us to contrast the estimation of marginal quantities
$\gamma_{0,M}$ and $\gamma_{1,M}$ and also those for subgroups of participants
based on baseline covariates $X$ from models 2 and 4. For our purposes, we
first compare the marginal estimates for $\gamma_{1}$ and $\gamma_{0}$ which
are directly calculated in settings 1 and 3 and are derived by marginalization
for models 2 and 4. The definition of the endpoints (either the difference
from baseline $D$ and corresponding $\gamma_{.D}$ or the original value $O$
and $\gamma_{.O}$) is important as $\gamma_{1,D}\neq\gamma_{1,O}$ when there
is an effect of $X$. Further, the validity of conditional independence varies:
when $T(0)\perp S(1)|T(1),X$, it is also true that $T(0)^{D}\perp
S(1)|T(1)^{D},X.$ However, this corresponding satisfaction of conditional
independence does not necessarily hold between $T(0)\perp S(1)|T(1)$ and
$T(0)^{D}\perp S(1)|T(1)^{D}$ when we do not condition on $X$. To quantify
this discrepancy when this condition is not met, we calculate the deviation
from meeting the conditional independence requirement, i.e. how incorrect it
is to make this conditional independence assumption (see Web Table 1 in the
supporting information).

###  4.1 Simulation Results

We present the posterior mean point estimates averaged over 100 datasets using
the observed data algorithm. To assess variability, we also report the
standard deviation of the Bayesian point estimates and the average standard
error of the parameter value (the posterior standard deviation) within each
replication. We run each MCMC for 3,000 iterations and ensure convergence is
reached using traceplots. Since our goal is to effectively validate surrogate
endpoints, we focus on inference for the quantities $\gamma_{0}$ and
$\gamma_{1}$. The results show the contrast in estimation accuracy and
efficiency when the constraints are enforced during estimation compared to
when no assumptions are imposed during the MCMC procedure.

###  4.2 Marginal Estimates

Below [Figure 1] (and in Web Table 1) are results with sample size $n=100$
using a truncated Beta($5,6$) prior. Bias for rows 1 and 2 as well as 3 and 4
is calculated based on the same true, marginal values. Overall, the identified
mean and variance parameters are estimated with little bias. However, the
nonidentified correlations are sensitive to the prior distribution. Compared
to Web Table 1, for comparability across settings Figure 1 has been adjusted
by the empirical variability in the validation estimates from complete
counterfactual data: we calculated the maximum likelihood estimates under a
scenario where we would observe all $n$ counterfactual outcomes and fit a
model for $T(1)-T(0)$ conditional on $S(1)$. We adjusted the bias and standard
errors relative to the standard deviation of the estimates across 100
simulations to provide quantities that can interpreted in proportion to the
amount of variability in the data and the estimates.

These results show there is often reduced bias in the estimates when making
conditional independence assumptions depending on the exact setting and prior
specification. Notably, the credible intervals for the nonidentified
parameters are very conservative for all settings as expected due to the non-
identified correlation parameter and its associated relatively weak prior, and
the corresponding coverage probabilities are near one. Further, the within
simulation average, over-estimation of the standard error decreases when
making conditional independence assumptions, though the SD of the estimates is
not necessarily smaller when making these assumptions (through less under-
identification in the conditional independence model, this reduction in SE
improves the agreement between the average SE and SD). We see that
conditioning on $X$ can also decrease the standard error of the estimates,
particularly in the scenario where we do not make conditional independence
assumptions.

###  4.3 Sensitivity to Distributional Assumptions

We show results in supporting information Web Table 2 for fitting conditional
models when the outcomes are heavy tailed (t-distributed) or skewed (gamma-
distributed). When data follows these non-Gaussian distributions, the
estimation results do vary from previous simulations. While there is some
increased bias and variability of the estimates as compared to those with a
normal distribution, generally there is some robustness in the estimation such
that the conclusions regarding surrogate validity seem to be similar in the
settings we considered.

###  4.4 Subgroup Analysis

Now we focus on the results of simulation setting $E$ where subgroups exist
for the models fit using methods 2 and 4 (conditional on $X$). Since $X$ is
binary in this setting, we report the conditional values of $\gamma_{0}$ when
$X=0$ and $X=1$. [Table 2] Whereas Figure 1 shows marginalized estimates of
the quantities, Table 2 shows the estimates of $\gamma_{0}$ and $\gamma_{1}$
conditional on the values of $X$, which allows us to determine if the
surrogate is valid for patients with certain baseline characteristics. Looking
at the validity of $S$ as a surrogate for subgroups in setting $E$, we see
that the surrogate is valid for $X=0$ but not for $X=1$. If $X$ were gender,
for example, this would indicate that the surrogate is valid only for males
and not females. However, this is only detected for the model where we make a
conditional independence assumption. When we do not, the credible interval of
$\gamma_{1}$ overlaps with 0 in almost all simulation replications, signaling
that we are not able to declare the surrogate as valid.

##  5 Duchenne Muscular Dystrophy Data Example

We generate “observed" data for three normally-distributed outcomes aimed to
mimic an ongoing clinical trial and the natural progression and deterioration
of ambulatory function for Duchenne muscular dystrophy patients.29 Since the
trial has not been unblinded, the parameter values were chosen to match the
preliminary data when available, and for other parameters the values were
chosen in a subjective way to be what the authors considered as reasonable.
Based on the literature and company analysis standards, to assess surrogacy,
we condition on age at baseline, $A$, and the baseline NSAA score measurement,
$X$.

###  5.1 Data Example Results

First we show the estimated treatment effects that would be calculated in a
clinical trial measuring the efficacy of the treatment based on the true,
observed endpoints $T$. [Table 3] It is clear that the standard error
estimates for $\gamma_{0}$ and $\gamma_{1}$ are markedly smaller, and credible
intervals more narrow, when we make the conditional independence assumptions.
We can see that marginally, the surrogate of micro-dystrophin is not a valid
surrogate for improvement of the NSAA score across the entire study sample.
However, there are strong effects of both age and baseline measurement on the
outcome. We can identify a region of the covariate space based on age and
baseline measurement where the surrogate is valid for a subgroup of patients.
We create CEP curves for these data to show what the surface looks like when
we stratify based on covariates. [Figure 2]

The CEP plot is conditional on baseline NSAA and age, and we demonstrate that
the surrogate is valid for those at four years of age. Since age is modeled
with linear and quadratic effects, the surrogate will be invalid for those six
and older as the estimated CEP curve moves farther from 0. This demonstrates
that due to the combination of natural growth and degeneration due to disease
over time, the surrogate would only be valid within a certain younger patient
population.

We also explore the consequences of different prior distributions for the non-
identified parameter $\theta_{T}$ and compare this to fixing $\theta_{T}$ at
some value, shown in Web Figure 1. We see that the results for $\gamma_{0}$
and $\gamma_{1}$ are somewhat sensitive to the choice, but the conclusion that
the surrogate is valid holds for values of $\theta_{T}$ which we believe to be
reasonable. The results are the same for both the imputation and observed data
algorithms.

##  6 Discussion

In this work, we have focused on incorporating baseline covariates into the
validation process for surrogate endpoints. Considering the harmful
implications of incorrectly validating a surrogate endpoint, it may be
worthwhile to consider the CEP as a function of $X$ and identify potential
subgroups of patients for which a surrogate is appropriate. It is important to
assess the plausibility of conditional independence assumptions even in this
context of adjusting for baseline covariates which we believe makes the
assumptions more likely. Introducing these conditional independence
constraints aids in improved estimation properties, but they should not be
implemented without proper, context-dependent reasoning. Further, there are
many ways to implement the constraints, although these simulations suggest
certain strategies (using observed data only) may reduce the burden of
imputing potential outcomes but still require reasonably well-specified prior
distributions.

There are many ways in which these methods and simulations will be extended,
particularly to verify its robustness. Work is ongoing to incorporate baseline
covariates and conditional independence assumptions into this framework while
fully utilizing the longitudinal data.

## Data Availability Statement

Data sharing is not applicable to this article as no new data were created or
analyzed in this study.

## Acknowledgements

This work was supported by the National Science Foundation (DGE 1256260) to
ER; the NIH (CA129102) to JMGT. Thanks to Walter Dempsey for useful
discussions.

## References

  1. 1.

Prentice R. Surrogate endpoints in clinical trials: definition and operational
criteria. Statistics in Medicine. 1989; 8(4), 431-40.

  2. 2.

VanderWeele TJ. Surrogate measures and consistent surrogates. Biometrics.
2013; 69(3), 561-565.

  3. 3.

Rubin DB. Estimating causal effects of treatments in randomized and
nonrandomized studies. Journal of Educational Psychology. 1974; 66(5),
688-701.

  4. 4.

Little R, Rubin D. Causal effects in clinical and epidemiological studies via
potential outcomes: concepts and analytical approaches.  Annual Review of
Public Health. 2000; 21(1), 121-145.

  5. 5.

Frangakis C, Rubin D. Principal stratification in causal inference.
Biometrics. 2002; 58(1), 21-29.

  6. 6.

Gilbert P, Hudgens M. Evaluating candidate principal surrogate endpoints.
Biometrics. 2008; 64(4), 1146-1154.

  7. 7.

Mendell JR, Sahenk Z, Lehman K, Nease C, Lowes LP, Miller NF, et al.
Assessment of Systemic Delivery of rAAVrh74. MHCK7. micro-dystrophin in
Children With Duchenne Muscular Dystrophy: A Nonrandomized Controlled Trial.
JAMA Neurology. 2020; 77(9), 1122-1131.

  8. 8.

Halloran ME, Préziosi MP, Chu H. Estimating vaccine efficacy from secondary
attack rates. Journal of the American Statistical Association. 2003; 98(461),
38-46.

  9. 9.

Préziosi MP, Halloran ME. Effects of pertussis vaccination on disease: vaccine
efficacy in reducing clinical severity. Clinical Infectious Diseases. 2003;
37(6), 772-779.

  10. 10.

Hudgens MG, Halloran ME. Causal vaccine effects on binary postinfection
outcomes.  Journal of the American Statistical Association. 2006; 101(473),
51-64.

  11. 11.

Follmann D. Augmented designs to assess immune response in vaccine trials.
Biometrics. 2006; 62(4), 1161-1169.

  12. 12.

Gabriel EE, Gilbert PB. Evaluating principal surrogate endpoints with time-to-
event data accounting for time-varying treatment efficacy.  Biostatistics.
2014; 15(2), 251-265.

  13. 13.

Gabriel EE, Follmann D. Augmented trial designs for evaluation of principal
surrogates.  Biostatistics. 2016; 17(3), 453-467.

  14. 14.

Gilbert PB, Qin L, Self SG. Evaluating a surrogate endpoint at three levels,
with application to vaccine development.  Statistics in Medicine. 2008; 62(4),
4758-4778.

  15. 15.

Wolfson J, Gilbert P. Statistical identifiability and the surrogate endpoint
problem, with application to vaccine trials. Biometrics. 2010; 66(4),
1153-1161.

  16. 16.

Huang Y, Gilbert PB, Wolfson J. Design and estimation for evaluating principal
surrogate markers in vaccine trials.  Biometrics. 2013; 69(2), 301-309.

  17. 17.

Zhuang Y, Huang Y, Gilbert PB. Simultaneous Inference of Treatment Effect
Modification by Intermediate Response Endpoint Principal Strata with
Application to Vaccine Trials. The International Journal of Biostatistics.
2019; 16(1).

  18. 18.

Gilbert PB, Huang Y. Predicting overall vaccine efficacy in a new setting by
re-calibrating baseline covariate and intermediate response endpoint effect
modifiers of type-specific vaccine efficacy.  Epidemiologic Methods. 2016;
5(1), 93-112.

  19. 19.

Zigler CM, Belin TR. A Bayesian approach to improved estimation of causal
effect predictiveness for a principal surrogate endpoint. Biometrics. 2012;
68(3), 922-932.

  20. 20.

Conlon AS, Taylor JM, Elliott MR. Surrogacy assessment using principal
stratification when surrogate and outcome measures are multivariate normal.
Biostatistics. 2014; 15(2), 266-283.

  21. 21.

Alonso A, Van der Elst W, Meyvisch P. Assessing a surrogate predictive value:
a causal inference approach. Statistics in medicine. 2017; 36(7), 1083-1098.

  22. 22.

Conlon A, Taylor J, Li Y, Diaz-Ordaz K, Elliott M. Links between causal
effects and causal association for surrogacy evaluation in a gaussian setting.
Statistics in Medicine. 2017; 36(27), 4243-4265.

  23. 23.

Taylor JM, Conlon AS, Elliott MR. Surrogacy assessment using principal
stratification with multivariate normal and Gaussian copula models. Clinical
Trials. 2015; 12(4), 317-322.

  24. 24.

Ding P, Feller A, Miratrix L. Randomization inference for treatment effect
variation.  arXiv preprint. 2014; arXiv:1412.5000.

  25. 25.

Parast L, Cai T, Tian L. Nonparametric estimation of the proportion of
treatment effect explained by a surrogate marker using censored data.
Technical Report. 2016.

  26. 26.

Parast L, McDermott MM, Tian L. Robust estimation of the proportion of
treatment effect explained by surrogate marker information. Statistics in
Medicine. 2016; 35(10) 1637-1653.

  27. 27.

Barnard J, McCulloch R, Meng XL. Modeling covariance matrices in terms of
standard deviations and correlations, with application to shrinkage.
Statistica Sinica. 2000; 1281-1311.

  28. 28.

Gustafson P. What are the limits of posterior distributions arising from
nonidentified models, and why should we care?. Journal of the American
Statistical Association. 2009; 104(488), 1682-1695.

  29. 29.

Muntoni F, Domingos J, Manzur AY, Mayhew A, Guglieri M, UK NorthStar Network,
et al. Categorising trajectories and individual item changes of the North Star
Ambulatory Assessment in patients with Duchenne muscular dystrophy. PloS One.
2019; 14(9), e0221097.

###  6.1 Tables and Figures

| A | B | C | D | E
---|---|---|---|---|---
| Valid $S$ | Valid $S$: No X Effect | Invalid $S$ | Valid $S$ for SG | Valid $S$ for SG
$\sigma_{X}$ | 0.5 | 0.5 | 0.5 | 0.5 | NA
$\delta_{4}$ | 1 | 1 | 1 | 1 | NA
$\omega_{1}$ | 2 | 2 | 2 | 2 | 2
$\omega_{2}$ | 0 | 0 | 0 | 0 | 0
$\omega_{3}$ | 3 | 3 | 3 | 3 | 3
$\omega_{4}$ | 1 | 0 | 1 | 3 | -0.75
$\omega_{5}$ | 4.1 | 4.1 | 4.1 | 4.1 | 4.1
$\omega_{6}$ | 1 | 0 | 1 | 1 | 2
$\epsilon_{S1}=\epsilon_{T0}=\epsilon_{T1}$ | 1 | 1 | 1 | 1 | 1
$\theta_{10}$ | 0.15 | 0.15 | 0.15 | 0.08 | 0.15
$\theta_{11}$ | 0.7 | 0.7 | 0.7 | 0.3 | 0.7
$\theta_{T}$ | 0.21 | 0.21 | 0.21 | 0.26 | 0.21
$\gamma_{0,O}$ | 0 | 0 | -1.00 | -1.35 | 1.31
$\gamma_{1,O}$ | 0.55 | 0.55 | 0.55 | 0.22 | 0.58
$\gamma_{0,D}$ | -0.06 | 0 | -1.02 | -1.33 | NA
$\gamma_{1,D}$ | 0.58 | 0.55 | 0.56 | 0.22 | NA

Table 1: Generative parameter values for the five scenarios to compare
definition of the endpoint and using baseline covariates. SG stands for
subgroup based on $X$. The meaning of these parameters is shown in equation 2.
Note ‘difference from baseline’ (subscript $D$) is not defined when $X$ is
binary and $T$ is continuous ($E$).

Setting | Fit Conditional | $\gamma_{0,C}$ |  |  |  |  | $\gamma_{0,C}$ |  |  |  |  | $\gamma_{1,C}$ |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| Independence | True | Est | SE | SD | Covers | True | Est | SE | SD | Covers | True | Est | SE | SD | Covers
| Assumption | $X=0$ |  |  |  | 0 | $X=1$ |  |  |  | 0 |  |  |  |  | 0
2E | $T(0)\perp S(1)|T(1),X$ | 0 | 0.064 | 0.472 | 0.308 | 1.00 | 2.75 | 2.778 | 0.470 | 0.312 | 0.00 | 0.55 | 0.512 | 0.179 | 0.089 | 0.00
2E | None | 0 | 0.078 | 1.024 | 0.309 | 1.00 | 2.75 | 2.784 | 1.013 | 0.312 | 0.04 | 0.55 | 0.509 | 0.481 | 0.089 | 0.98

Table 2: Simulation results demonstrating effect of estimating subgroups.
Estimates of $\gamma_{0}$ and $\gamma_{1}$ are conditional on the values of
$X$, so values of $\gamma_{0}$ estimates are conditional on $X=0,1$. The
columns denoting ‘covers 0’ indicate what proportions of simulations have
credible intervals that do contain 0. This helps determine for how many
simulations the surrogate would be considered valid. For $\gamma_{0}$, the
credible interval covering 0 denotes a valid surrogate, while the interval of
$\gamma_{1}$ should not. Figure 1: Simulation results for the two definitions
of the endpoint (here categorized by superscript to denote $O$ = Original
Endpoint (Settings 1 and 2) compared to $D$ = Difference from baseline
(Settings 3 and 4)) and different generating parameter values $(A-E)$. The
values shown below are the bias and variability (both the average within-
sample standard error across datasets and the standard deviation of the point
estimates) of the validation quantities that are adjusted by the variability
of the hypothetical surrogate values if all counterfactuals were to be
observed (from the full data).

Endpoint | Treatment | Conditional Independence | $\gamma_{0}$ | $\gamma_{0}$ | $\gamma_{1}$ | $\gamma_{1}$
---|---|---|---|---|---|---
Type | Effect | Assumption | Estimate | SE | Estimate | SE
$T^{D}$ | 0.227 | $T^{D}(0)\perp S(1)|T^{D}(1),X,A$ | -3.040 | 0.393 | 0.393 | 0.101
|  | None | -2.980 | 2.235 | 0.386 | 0.279
$T$ | 0.271 | $T(0)\perp S(1)|T(1),X,A$ | -3.046 | 0.826 | 0.394 | 0.102
|  | None | -2.986 | 2.232 | 0.386 | 0.279

Table 3: Simulated Muscular dystrophy estimates of treatment effect and
marginal surrogacy quantities. Figure 2: The CEP plot shows the conditional
functions of $\gamma_{0}$ and $\gamma_{1}$ across possible values of $S(1)$.
The empirical distribution of $S(1)$ is shown in the blue density curve.

.