# THE FIRST COMPREHENSIVE DATASET WITH MULTIPLE DISTORTION TYPES FOR VISUAL
JUST-NOTICEABLE DIFFERENCES

###### Abstract

Recently, with the development of deep learning, a number of Just Noticeable
Difference (JND) datasets have been built for JND modeling. However, all the
existing JND datasets only label the JND points based on the level of
compression distortion. JND models learned from such datasets can only be used
for image/video compression. Hence, a generalized JND modeling should take
more kinds of distortion types into account. To benefit JND modeling, this
work establishes a generalized JND dataset with a coarse-to-fine JND
selection, which contains 106 source images and 1,642 JND maps, covering 25
distortion types. To this end, we proposed a coarse JND candidate selection
scheme to select the distorted images from the existing Image Quality
Assessment (IQA) datasets as JND candidates instead of generating JND maps
ourselves. Then, a fine JND selection is carried out on the JND candidates
with a crowdsourced subjective assessment.

Index Terms—  Just Noticeable Difference, Human Visual System, Mean Opinion
Scores (MOS), perception modeling, dataset

##  1 Introduction

JND is a metric for assessing the visual redundancy of the HVS, and it has
been widely applied to computer vision and multimedia signal processing
applications such as perceptual image and video compression [1] [2], image
enhancement [3], watermarking [4], and so on. JND modeling, which has been
studied for many years, tries to precisely predict the visual redundancy of
the HVS for a given visual content. Traditional JND models [5] [6] tried to
predict the JND threshold for each pixel or each coefficient of the sub-bands
based on the features of the HVS and their associated maskings.

Fig. 1: Illustration of some source images in our dataset

Recently, a few works tried to extend the JND concept to a picture level and
proposed the picture-wise JND. To utilize the powerful deep learning
techniques in the JND modeling, lots of JND datasets [7, 8, 9] were built. Jin
et al. [7] created the MCL-JCI dataset, which consists of 5,000 corresponding
JPEG distorted versions with quality factors (QF) ranging from 1 to 100 and 50
source images in the $1920\times 1080$ resolution. They conducted subjective
quality assessment tests involving 150 participants watching the source image
and a compressed image side by side on a TV. JND samples were gathered from 30
individuals for a specific image. Then, Shen et al. [8] used 202 source images
and 7,878 encoded versions of those images with a resolution of $1920\times
1080$ to create a JND dataset based on the impending video coding standard
Versatile Video Coding (VVC). Each source image was compressed by VTM5.0
coding with quantization parameters (QP) ranging from 13 to 51. The subjective
tests were conducted in a carefully monitored lab setting. 20 subjects
evaluated 20 PJND samples for each source image. Lin et al. [9] expanded the
datasets mentioned above and built a large PJND dataset called KonJND-1k. They
conducted subjective JND assessment studies using the flicker test via
crowdsourcing instead of in the laboratory environment. This dataset contains
1,008 source images as well as distorted versions created using JPEG and BPG
compression. The study involved 503 workers in total, resulting in 61,030 PJND
samples and an average of 42 samples per source image. However, all of these
JND datasets in existence only took into consideration the impacts of
compression distortion without considering the effects of the other distortion
types. That is, the JND points were labeled based on the levels of compression
distortion, _e.g._ , Quality Factor (QF), Quantization Parameter (QP), _etc_.
Therefore, JND models learned from such kinds of datasets were limited to
compression-relevant applications. In order to predict a JND model so that it
has more wide applications, establishing a generalized JND dataset that covers
all existing distortion types is necessary and significant. However, the
challenging part is how to obtain distorted images that cover all distortion
types with different distortion levels so that JND maps can be selected from
them.

Recently, many IQA datasets [10, 11, 12, 13] have been established, which
contain source images and their associated distortion images as well as
subjective scores, _e.g._ , Mean Opinion Scores (MOS). Commonly, their
distortion images cover various distortion types under various distortion
levels. In view of this, we can directly pick out the JND maps from the
distorted images in the IQA datasets. On the one hand, the generation of
distortion images that cover all distortion types with different distortion
levels can be saved. On the other hand, as the MOS reflects the perceptual
quality of distorted images, it can be used for selecting JND maps. Therefore,
we directly select the JND maps from IQA datasets with a proposed coarse-to-
fine JND map selection scheme, which contains two main steps: 1) coarse JND
candidate selection with MOS as the threshold and 2) fine JND map selection
with crowdsourced subjective assessment. The first step helps us make a fast
JND candidate selection from huge distorted images, which can save time and
human resources in subsequent subjective assessment. The second step
subjectively confirms the final JND maps from the selected JND candidates in
the first step by the subjects. Finally, we establish a generalized JND
dataset, containing 106 source images (part of them are shown in Fig. 1) and
1,642 JND maps, with 25 distortion types being involved. It should be
mentioned that this is the first JND dataset that covers various distortion
types.

##  2 Coarse-to-Fine JND Selection

###  2.1 Coarse JND Candidate Selection

To select JND images from these datasets, we use the MOS of distorted images
as the preliminary selection criterion, since the MOS is a subjective viewing
score, reflecting the perceptual quality of distorted images among all the
subjects. Commonly, the higher the distortion image quality, the smaller the
difference between it and its associated source image, and vice versa.
Therefore, we can use MOS to select distorted images with high perceptual
quality as the JND maps, so as to reduce the amount of data in subsequent
subjective tests. The involved IQA datasets in this work include TID2013 [13],
TID2008 [12], and KADID-10k [10], which have a total of 14,825 distortion
images from 106 source images, and 39 distortion types, each with 4 or 5
distortion levels.

Fig. 2: Histogram of the distribution of normalized MOS in TID2008, TID2013,
and KADID-10k datasets, which shows that different IQA datasets have different
MOS distributions.

Since different IQA datasets use different algorithms to subjectively evaluate
distortion images, the range of the obtained MOS values is different. In view
of this, we first normalize the MOS values. This process can be formulated as

$Q=\frac{P-Min}{Max-Min},$ (1)

where $P$ represents the MOS value of the distortion image in a certain IQA
dataset, $Min$ and $Max$ are the minimum and maximum values of the MOS in the
current IQA dataset. $Q$ represents the normalized MOS. The range of the
normalized MOS is [0,1]. Fig. 2 shows the MOS distribution of each IQA dataset
after normalization, and it can be seen that the MOS distribution in different
datasets is different. Among them, the MOS distribution of TID2013 and TID2008
conforms to the normal distribution, while the KADID-10k dataset is relatively
evenly distributed in each interval. To determine a threshold that can
distinguish between JND and non-JND maps, we randomly sample 10%-20% of the
source images and all their corresponding distorted images from each IQA
dataset and then pick out JND maps from them through a subjective viewing test
by well-trained and experienced subjects. We only sample 10%-20% of the whole
images because the number of images is huge and we don’t need to ensure the
high accuracy of the threshold due to the fine JND map selection after. As a
result, 112 JND maps are confirmed out of 2,190 distorted images with MOS
above 0.7. To determine the MOS threshold, we calculate the average MOS of the
confirmed JND images for each image dataset and then reduce the calculated
average MOS by about 5% so that more JND-eligible maps can be involved. Then,
we obtain the MOS thresholds 0.855 for the KADID-10k dataset, and 0.795 for
the TID2013 and TID2008 datasets, respectively. Afterward, we use the MOS
thresholds above to pick out images from their corresponding datasets as the
JND candidates. In this way, we coarsely obtain 1,962 JND candidates.

###  2.2 Fine JND Map Selection

After coarse JND selection, there are still some JND candidates that can be
noticed the differences compared with their associated source images. In view
of this, we conduct a crowdsourced subjective assessment of the JND candidates
to achieve fine JND selection by abandoning non-JND maps. The crowdsourced
subjective assessment was conducted on Amazon Mechanical Turk (AMT). On the
AMT platform, requesters—companies, organizations, or individuals—create and
submit tasks requiring human intelligence for workers, who may be paid for
each task successfully completed.

In this work, we adopt the flicker test as used in [14], where the JND
candidate and its associated source image are toggled back and forth
successively at a frequency of 8 Hz. Here, we set three options for subjects:
obvious flicker, slight flicker, and no flicker. The definitions of the three
options are as follows. The obvious flicker indicates that the flicker is
intense, affecting the entire image. The slight flicker means that slight
changes can be seen in the image, such as a subtle variation in the image
brightness. These changes can resemble a mild film grain in a static scene of
a movie and are often localized to small regions of the image. No flicker
represents the image is static, there is no visible change in the image, no
matter how small it is.

According to the definition of JND as aforementioned in Sec. 1, no flicker and
slight flicker JND candidates are selected as the final JND maps. Here, we
randomly divide all the 1962 JND candidates into 103 groups. Each group
contains 19 or 20 JND candidates, which are assessed by 30 workers (subjects).
A total of $1962\times 30$ results are obtained.

Considering that the screens used by the workers may have different sizes and
resolutions, resulting in different physical sizes of pictures displayed on
the screen, we carry out the calibration process in [9] to make the image
appear the same physical size on all workers’ computers. That is, workers must
prepare a card the size of a credit card ($85mm\times 53.98mm$) and change the
size of a frame on the screen to fit the card. By computing the Logical Pixel
Density (LPD) of the display in Pixels Per Inch (PPI), we can estimate the
physical size of the display. For more details refer to [9]. Following the
calibration, we instruct the workers to change their viewing distance to
$30cm$ in accordance with the trigonometric calculation [15], [16] and ISO
standard [17].

After passing a training phase, the workers will enter the assessment phase.
In the assessment phase, 19 or 20 images are set for each task, and the
workers choose the most suitable option from the three options of no flicker,
slight flicker, and obvious flicker, as aforementioned.

##  3 Results Processing and Analysis

###  3.1 Outlier Removal

To get effective results, we need to remove outliers in the collected results.
Since we can not guarantee that all workers complete the test with
seriousness, for example, some workers may randomly choose an option during
the test, or workers may submit unreliable answers because they have been
tested for a long time. In this paper, we adopt the result processing method
mentioned in the ITU-R Recommendation BT.500 [18] to remove unreliable
results. We use the scores 1, 2, and 3 to represent no flicker, slight
flicker, and obvious flicker respectively. Then, we calculate the mean and
standard deviation of the scores rated by all the subjects for each image as
follows

$\bar{u}_{j}=\frac{1}{N}\sum_{i=1}^{N}{u_{ij}},$ (2)
$S_{j}=\sqrt{\sum_{i=1}^{N}{\frac{\left(\bar{u}_{j}-u_{ij}\right)^{2}}{\left(N-1\right)}}}.$
(3)

$u_{ij}$ denotes the score rated by the $i$-th subject on the $j$-th image,
where $i=1,2,...,N$ and $j=1,2,...,M$. $\bar{u}_{j}$ and $S_{j}$ are the mean
and standard deviation of scores rated by all the subjects for the $j$-th
image. Subsequently, we adopt the confidence interval in [18] for outlier
removal, that is

$\left[\bar{u}_{j}-\delta_{j},\bar{u}_{j}+\delta_{j}\right],\text{where}\
\delta_{j}=1.96\frac{S_{j}}{\sqrt{N}}.$ (4)

$[\cdot]$ is rounding operation. $u_{ij}$ outside such an interval is
considered as an outlier. After removing the outliers, we re-calculate the
average score of each image. We regard the image with an average score of less
than 2.5 as a JND map.

###  3.2 Result Analysis

Table 1: Numbers of JND Maps and The Ratio of The Number of JND Maps to The Total Samples for Different Distortion Types No. Distortion types | Num. | Ratio
---|---|---
01 Additive Gaussian noise | 19 | 0.030
02 Additive noise in color components | 115 | 0.183
03 Masked noise | 50 | 0.222
04 High frequency noise | 38 | 0.167
05 Denoise | 50 | 0.079
06 Multiplicative noise | 55 | 0.104
07 Comfort noise | 18 | 0.144
08 Motion blur | 91 | 0.225
09 Gaussian blur | 128 | 0.203
10 Color diffusion | 60 | 0.148
11 Color shift | 25 | 0.062
12 Color quantization | 25 | 0.047
13 Color saturation | 70 | 0.132
14 JPEG | 95 | 0.151
15 JPEG transmission errors | 26 | 0.116
16 JPEG2000 | 123 | 0.195
17 Darken | 138 | 0.341
18 Brighten | 68 | 0.168
19 Mean shift | 233 | 0.370
20 Jitter | 40 | 0.099
21 Non-eccentricity patch | 25 | 0.040
22 Quantization | 36 | 0.089
23 Contrast change | 54 | 0.086
24 Sparse sampling and reconstruction | 13 | 0.104
25 Chromatic aberrations | 47 | 0.376

As a result, we select a total of 1,642 distortion images (corresponding to
106 source images) as the JND maps, which contain 25 distortion types. Table 1
shows the 25 distortion types, their number of JND maps, and the ratio of the
number of JND maps to the total number of samples. From the table, we can see
that the number of JND maps with different distortion types varies greatly.
There are two reasons to explain this. First, since we select JND maps from
three IQA datasets KADID-10k, TID2008 and TID2013, some distortion types are
not contained by all three datasets such as Comfort noise and Sparse sampling
and reconstruction, which causes different numbers of samples for different
distortion types. The distortion type with a small number of samples probably
has fewer JND maps and vice versa. Second, the HVS has different sensibilities
to different types of distortion [6]. For the distortion type with a high
sensibility of the HVS, even low-level distortion may be perceived by the HVS
and this will make fewer JND maps can be obtained from the IQA datasets, and
vice versa. Therefore, the distortion with a higher sensibility to the HVS,
the fewer JND maps we obtained from the samples with such type of distortion.
Besides, for each distortion type, we also calculate the ratio of the number
of JND maps and total samples. The difference in such ratios reflects a
certain extent to which types of distortion are more easily sensitive to the
HVS.

A good JND dataset tends to contain a wide variety of visual content. To
demonstrate the diversity of the JND maps in our dataset, we calculated the
spatial information (SI) [19] and the colorfulness (CF) [20] of 106 source
images, as shown in Fig. 3. The spatial information reflects the spatial
complexity of the image, and the colorfulness reflects the richness of the
color of the image. In Fig. 3, we can see that the selected source images do
have wide coverage in this plot. Furthermore, we divide the 106 source images
into 8 semantic categories, as shown in Table 2.

Fig. 3: Distribution of spatial information and colorfulness of source images in our dataset. Table 2: Number of Source Images in Each Semantic Category Content | Number
---|---
People | 12
Animals | 16
Plants | 12
Objects | 21
Landscape | 9
Outdoor | 15
Building | 12
Nature | 9

##  4 Conclusions

A JND dataset containing 106 source images and 1,642 JND maps with 25
distortion types involved has been established in this work. To have such a
dataset, we coarsely select the JND candidates from IQA datasets through the
MOS thresholds, and then subjectively select the final JND map from the JND
candidates by conducting crowdsourced subjective assessments, which can
greatly save time and human resources spent generating distorted images by
ourselves. This is the first JND dataset covering multiple distortion types,
which can be used to learn JND models for more distortion types, not just
limited to compression distortion types contained in existing JND datasets,
and thus have a wider range of applications. For example, we can use it to
build a JND model to optimize the image transmission process and reduce the
image distortion caused by packet loss during data transmission. Furthermore,
considering that pictures on social media now contain multiple distortions,
this dataset can also be used to build a more generalized JND model, which can
well simulate the perception of human vision for any type of distortion, thus
revealing the characteristics of the human visual system and improving the
quality of images on social media sites. Besides, our dataset has strong
expansibility. In our future work, we will extend this dataset by applying the
proposed coarse-to-fine JND selection scheme to more IQA datasets such as
PIPAL [21], CSIQ [11], and KADIS-700k [22].

## References

  * [1] H. R. Wu, A. R. Reibman, W. Lin, F. Pereira, and S. S. Hemami, “Perceptual visual signal compression and transmission,” _Proceedings of the IEEE_ , vol. 101, no. 9, pp. 2025–2043, 2013.
  * [2] J. Kim, S.-H. Bae, and M. Kim, “An hevc-compliant perceptual video coding scheme based on jnd models for variable block-sized transform kernels,” _IEEE Transactions on Circuits and Systems for Video Technology_ , vol. 25, no. 11, pp. 1786–1800, 2015.
  * [3] Z. Cheng, H. Sun, M. Takeuchi, and J. Katto, “Performance comparison of convolutional autoencoders, generative adversarial networks and super-resolution for image compression,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , 2018, pp. 2613–2616.
  * [4] J. Li, H. Zhang, J. Wang, Y. Xiao, and W. Wan, “Orientation-aware saliency guided jnd model for robust image watermarking,” _IEEE Access_ , vol. 7, pp. 41 261–41 272, 2019.
  * [5] C.-H. Chou and Y.-C. Li, “A perceptually tuned subband image coder based on the measure of just-noticeable-distortion profile,” _IEEE Transactions on circuits and systems for video technology_ , vol. 5, no. 6, pp. 467–476, 1995\.
  * [6] J. Wu, L. Li, W. Dong, G. Shi, W. Lin, and C.-C. J. Kuo, “Enhanced just noticeable difference model for images with pattern complexity,” _IEEE Transactions on Image Processing_ , vol. 26, no. 6, pp. 2682–2693, 2017.
  * [7] L. Jin, J. Y. Lin, S. Hu, H. Wang, P. Wang, I. Katsavounidis, A. Aaron, and C.-C. J. Kuo, “Statistical study on perceived jpeg image quality via mcl-jci dataset construction and analysis,” _Electronic Imaging_ , vol. 2016, no. 13, pp. 1–9, 2016.
  * [8] X. Shen, Z. Ni, W. Yang, X. Zhang, S. Wang, and S. Kwong, “A jnd dataset based on vvc compressed images,” in _2020 IEEE International Conference on Multimedia & Expo Workshops (ICMEW)_, 2020, pp. 1–6.
  * [9] H. Lin, G. Chen, M. Jenadeleh, V. Hosu, U.-D. Reips, R. Hamzaoui, and D. Saupe, “Large-scale crowdsourced subjective assessment of picturewise just noticeable difference,” _IEEE transactions on circuits and systems for video technology_ , vol. 32, no. 9, pp. 5859–5873, 2022.
  * [10] H. Lin, V. Hosu, and D. Saupe, “Kadid-10k: A large-scale artificially distorted iqa database,” in _2019 Eleventh International Conference on Quality of Multimedia Experience (QoMEX)_ , 2019, pp. 1–3.
  * [11] E. C. Larson and D. M. Chandler, “Most apparent distortion: full-reference image quality assessment and the role of strategy,” _Journal of electronic imaging_ , vol. 19, no. 1, pp. 011 006–011 006, 2010.
  * [12] N. Ponomarenko, V. Lukin, A. Zelensky, K. Egiazarian, M. Carli, and F. Battisti, “Tid2008-a database for evaluation of full-reference visual quality assessment metrics,” _Advances of Modern Radioelectronics_ , vol. 10, no. 4, pp. 30–45, 2009.
  * [13] N. Ponomarenko, L. Jin, O. Ieremeiev, V. Lukin, K. Egiazarian, J. Astola, B. Vozel, K. Chehdi, M. Carli, F. Battisti _et al._ , “Image database tid2013: Peculiarities, results and perspectives,” _Signal processing: Image communication_ , vol. 30, pp. 57–77, 2015.
  * [14] I. 29170-2, “Information technology – advanced image coding and evaluation – part 2: Evaluation procedure for visually lossless coding,” _British Standatds Institution_ , 2015.
  * [15] K. Rohrschneider, “Determination of the location of the fovea on the fundus,” _Investigative ophthalmology & visual science_, vol. 45, no. 9, pp. 3257–3258, 2004.
  * [16] B. V. Ehinger, K. Häusser, J. P. Ossandon, and P. König, “Humans treat unreliable filled-in percepts as more real than veridical ones,” _Elife_ , vol. 6, p. e21761, 2017.
  * [17] B. ISO and B. STANDARD, “Ergonomics of human-system interaction,” _British Standatds Institution_ , 2010.
  * [18] B. Series, “Methodology for the subjective assessment of the quality of television pictures,” _Recommendation ITU-R BT_ , pp. 500–13, 2012.
  * [19] T. Installations and L. Line, “Subjective video quality assessment methods for multimedia applications,” _Networks_ , vol. 910, no. 37, p. 5, 1999.
  * [20] S. Winkler, “Analysis of public image and video databases for quality assessment,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 6, no. 6, pp. 616–625, 2012.
  * [21] G. Jinjin, C. Haoming, C. Haoyu, Y. Xiaoxing, J. S. Ren, and D. Chao, “Pipal: a large-scale image quality assessment dataset for perceptual image restoration,” in _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XI 16_ , 2020, pp. 633–651.
  * [22] H. Lin, V. Hosu, and D. Saupe, “Deepfl-iqa: Weak supervision for deep iqa feature learning,” _arXiv preprint arXiv:2001.08113_ , 2020.


# Constructing a T-test for Value Function Comparison of Individualized
Treatment Regimes in the Presence of Multiple Imputation for Missing Data

Minxin Lu
Department of Biostatistics
University of North Carolina at Chapel Hill
North Carolina, U.S.A.
<EMAIL_ADDRESS>
&Annie Green Howard
Department of Biostatistics, Carolina Population Center
University of North Carolina at Chapel Hill
North Carolina, U.S.A.
<EMAIL_ADDRESS>
Penny Gordon-Larsen
Department of Nutrition, Carolina Population Center
University of North Carolina at Chapel Hill
North Carolina, U.S.A.
<EMAIL_ADDRESS>
&Katie A. Meyer
Department of Nutrition
University of North Carolina at Chapel Hill
North Carolina, U.S.A.
<EMAIL_ADDRESS>
&Hsiao-Chuan Tien
Carolina Population Center
University of North Carolina at Chapel Hill
North Carolina, U.S.A.
<EMAIL_ADDRESS>
&Shufa Du
Department of Nutrition, Carolina Population Center
University of North Carolina at Chapel Hill
North Carolina, U.S.A.
<EMAIL_ADDRESS>
&Huijun Wang
National Institute for Nutrition and Health
Chinese Center for Disease Control and Prevention
Beijing, China
<EMAIL_ADDRESS>
&Bing Zhang
National Institute for Nutrition and Health
Chinese Center for Disease Control and Prevention
Beijing, China
<EMAIL_ADDRESS>
&Michael R. Kosorok
Department of Biostatistics
University of North Carolina at Chapel Hill
North Carolina, U.S.A.
<EMAIL_ADDRESS>
###### Abstract

Optimal individualized treatment decision-making has improved health outcomes
in recent years. The value function is commonly used to evaluate the goodness
of an individualized treatment decision rule. Despite recent advances,
comparing value functions between different treatment decision rules or
constructing confidence intervals around value functions remains difficult. We
propose a t-test based method applied to a test set that generates valid
p-values to compare value functions between a given pair of treatment decision
rules when some of the data are missing. We demonstrate the ease in use of
this method and evaluate its performance via simulation studies and apply it
to the China Health and Nutrition Survey data.

_Keywords_ Value Function, T-test, Precision Medicine, Individualized
Treatment Regimes, Imputation

##  1 Introduction

The response to a particular treatment can vary among individuals due to the
influence of their unique characteristics. This variability in treatment
effect is particularly relevant in healthcare, where individualized treatment
regimes (ITRs) are employed. ITRs tailor treatments and preventive approaches
for individuals based on factors such as their socioeconomic status,
environment, lifestyle choices, and medical conditions. The application of
ITRs extends beyond healthcare and medical treatment assignments, encompassing
precision nutrition and behavioral interventions as well.

In these applications, decision-makers may seek to compare the effectiveness
of an ITR against, for example, the observed treatment assignment, a one-size-
fits-all treatment approach, or another ITR. This comparison is valuable for
both patients and healthcare providers in directing their efforts toward the
most effective treatment option. For example, individuals with or at risk for
hypertension are generally advised to increase their physical activity. The
American Heart Association, based on the Physical Activity Guidelines for
Americans 1, recommends that adults engage in 2.5-5 hours of moderate-
intensity physical activity per week; although it is noted that additional
health benefits can be achieved by achieving beyond 5 hours per week 2.
Nevertheless, this physical activity recommendation can be challenging to meet
for many people. In addition, it is possible that this level of physical
activity may not be beneficial for all individuals 3.

Thus, it is important to examine whether physical activity, especially at high
recommended levels, which might be challenging to implement, uniformly
benefits hypertension prevention and management across all demographic groups.
This may be particularly relevant in countries in which the prevalence of
anti-hypertensive medications is low, such as China. We aim to use ITRs to
differentiate between subpopulations that experience incremental benefits from
increasing their weekly physical activity to more (versus less or equal) than
5 hours per week and directly compare the ITR generated to population-level
recommendations. For subpopulations who are unable to adhere to recommended
activity levels, there are potential significant benefits conferred by
focusing on alternative behavior interventions, such as dietary modification.
Moreover, there are various methods for estimating the ITRs, such as
Q-learning 4 and D-learning 5. Direct comparisons of these ITRs are valuable
for determining whether one ITR results in a superior population-level
improvement, as compared to another ITR, or to a non-individualized treatment.

Observational data is frequently encountered in research studies and serves as
a vital data source for deriving ITRs. Observational epidemiologic studies
have several potential strengths, including increased representation of
variability in the target population, through large sample sizes and sampling
strategies to increase population representation. However, observational data
generally lack balance in treatment assignments. To address this confounding,
the average outcome for the population, assuming adherence to a specific
treatment regime, is weighted by inverse individual propensity scores, and
this weighted outcome is defined as the value function. Researchers have
commonly employed the value function to evaluate the performance of an ITR
derived from observational cohort data. The presence of missing data
constitutes a common challenge encountered within randomized trial datasets as
well as observational datasets. While large observation data potentially
provides more heterogeneity than smaller, more controlled studies,
observational studies are often more susceptible to missing data problems. One
efficacious approach involves the application of multiple imputation 6 whereby
multiple plausible imputed datasets are generated. Each imputed dataset
entails the substitution of missing values with imputed values derived from a
predictive model. Subsequently, standard statistical methods are applied to
each of the imputed datasets separately. This process yields a collection of
results, and the observed variability among these results across the imputed
datasets reflects the inherent uncertainty associated with missing values.
Multiple imputation allows comprehensive assessment of the robustness of the
study findings in the presence of missing data and has the potential to
enhance reliability and generalizability of the research outcomes.

Various approaches have been employed to estimate the variance of a value
function for ITRs, such as jackknife 7, cross-validation 8, 9, bootstrapping
10, and Q-function model-based approaches 11. Each is explained below.
However, these approaches tend to become intricate when using multiple
imputation and are not intended for directly comparing the value functions of
a pair of ITRs, such as estimating the variance around the difference between
two value functions. 7 employed the jackknife or leave-one-out cross-
validation method to estimate the variance of the value function and compared
24 individualized treatment regimes derived from 24 machine learning models.
This jackknife approach consistently estimates the variance around a value
function estimate and requires only weak assumptions such as requiring the
samples to be independent and identically distributed, and that as the sample
size becomes larger, the decision rule estimated from n-1 individuals
converges to the decision rule estimated from n individuals. However, it is a
time-consuming implementation and better suited for small datasets. Cross-
validation provides an alternative means to estimate variance. 8 and 9 used
cross-validation in their paper. The data is divided into $K$ folds, and the
ITR was estimated based on the $K-1$ folds and the value function was
evaluated based on the remaining fold. The process is repeated $K$ times to
get $K$ estimates of value functions and then the variance of value function
is calculated. A common choice of $K$ is $K=10$ which is called 10-fold cross-
validation. This method is less time-consuming than the jackknife, though it
has a tendency to understate the actual variance. This is because each data
point is utilized in both the training and testing sets, leading to a
correlation among the accuracy measures of each fold 12. Bootstrapping
provides an alternative approach to compute the variance of the value function
for an ITR. The standard n-out-of-n bootstrap method 10 involves randomly
selecting n observations with replacements from the original data with n
observations to create new datasets, which are then repeatedly generated 500
or 1000 times. Then, methods for finding the value function of the ITR are
applied to each sequence, resulting in slightly different value functions.
Subsequently, the variance and confidence interval of the ITR’s value function
are developed. Variations of the standard n-out-of-n bootstrap method include
double bootstrap 13, adaptive bootstrap 14, and m-out-of-n bootstrap 15. The
bootstrapping method requires repeating operations on the data, usually 500
times or more, making the approach computationally intensive. 11 have
developed a method for constructing the statistical inference of a policy’s
value function for reinforcement learning when either the number of decision
points or the sample trajectories diverge to infinity. Their approach involves
utilizing the sieve method to approximate the Q- function and employs
“SequentiAl Value Evaluation” to split data and iteratively find the optimal
policy and the value function estimate. Their method constructs a valid
confidence interval around the value function estimate that achieves the
nominal coverage. However, it requires estimation of the Q-function and relies
on the correctness of the Q-function model.

In this paper, we present a new method that enables the direct comparison of
any two ITRs via the value function and also provides a t-test-based p-value
for the significance of observed differences. Our approach addresses the
shortcomings described above in that: 1) it is less time-consuming to
implement than the jackknife or bootstrapping; 2) it circumvents the problem
of correlation in each fold experienced by cross-validation; and 3) it does
not rely on the estimation of the Q-function. Our method is suitable for both
observational studies and clinical trials. Moreover, our method enables
deriving variance estimates from multiple imputed datasets without the need
for additional replication or bootstrapping, making it particularly
advantageous and efficient when dealing with missing data. Specifically, our
approach provides a valid estimate of the variability surrounding both the
value function itself and the difference between the two value functions. With
our approach, it is possible to assess whether the optimal ITR significantly
outperforms the one-size-fits-all approach or another ITR estimated using a
different method. Furthermore, our method is characterized by its ease of
understanding and implementation, while providing theoretical guarantees for
estimator consistency and inference validity, including variance and
confidence interval calculations. For illustrative purposes, we use two
models, Q-learning 4 and D-learning 5 for ITR comparison, although more
sophisticated models can also be incorporated.

The following outlines the structure of this paper. Section 2 presents the
introduction of our method. In section 3, we showcase the application of our
method through simulation studies. Section 4 delves into the implementation of
our method on a real-world dataset. Finally, in Section 5, we discuss the
advantages and limitations of our approach.

##  2 Methods

We assume that the dataset is divided into a training set that estimates the
decision rule and a test set that evaluates this rule on new data, with $n$
and $m$ individuals respectively. The training data is independent of the test
data. The individual index is denoted by $i$, where $X_{i}$ represents the
covariate vector, $A_{i}$ the binary treatment assignment, and $Y_{i}$ the
outcome. Let $\hat{d}_{1,n}$ and $\hat{d}_{2,n}$ be two decision rules we want
to compare. Each $\hat{d}_{n}$ is a map from covariates to treatments. The
$\hat{d}_{1,n}$ and $\hat{d}_{2,n}$ can be zero-ordered decision rules (e.g.
assign everyone the same treatment), or ITRs as a function of the covariates
based on models from the training set (e.g. assign everyone who is younger
than 40 years old one treatment, and assign placebo otherwise). The variable
$n$ indicates that the rules are estimated using training data, and in case of
missing data, decision rules are derived from single or multiple imputation of
the training data. Our objective is to estimate:

  * •

The value function of the decision rule $d_{j,n}$ on the test set
$V_{m}(d_{j,n})=E(Y_{m}|A=\hat{d}_{j,n}(x))$, for j =1,2

  * •

The variance associated with the value function $Var(V_{m}(d_{j,n}))$

  * •

The difference between value functions of two decision rules
$V_{m}(d_{1,n})-V_{m}(d_{2,n})$

  * •

The variance of the difference between the value functions of two decision
rules $Var(V_{m}(d_{1,n})-V_{m}(d_{2,n}))$

  * •

The p-value of the t-test for the significance of the difference between the
value functions of two decision rules $V_{m}(d_{1,n})-V_{m}(d_{2,n})$

###  2.1 Propensity Score

Although randomized trial data are ideal for estimating the optimal ITR, many
researchers only have access to observational data. When using observational
data, the confounding effect is a main concern and can be mitigated by
propensity score modeling 16. The propensity score, denoted as
$\pi(A|X)=P(A|X)$, models the probability of receiving treatment $A$ given
covariate vector $X$. By weighting each individual with the inverse propensity
score, the observational data can approximate the data from a randomized
trial. Let $\hat{\pi}_{n}(a|x)$ be the estimated propensity score from the
training set, and assume $\hat{\pi}_{n}(a|x)=\pi(a|x,\hat{\theta}_{n})$, where
$\hat{\theta}_{n}$ are parameters for the propensity score model estimated
based on the training set. Let $\theta_{0}$ represent the true parameter
vector.

######  Assumption 1.

$\sqrt{n}(\pi(a|x,\hat{\theta}_{n})-\pi(a|x,\theta_{0}))=\sqrt{n}(\hat{\theta}_{n}-\theta_{0})^{T}\phi_{0}(a,x)+o_{p}(1),$
(1)

where $\sqrt{n}(\hat{\theta}_{n}-\theta_{0})\xrightarrow{N}(0,\Sigma_{0})$,
and $\Sigma_{0}$ is the limiting variance of
$\sqrt{n}(\hat{\theta}_{n}-\theta_{0})$. Let
$\Sigma_{n}=Var(\hat{\theta}_{n})$, then $\Sigma_{n}=\Sigma_{0}/n+o(1/n)$. And
$n\hat{\Sigma}_{n}$ is an estimate of $\Sigma_{0}$ obtained from the training
data.

######  Assumption 2.

Let $\hat{\phi}_{n}$ be an estimate of $\phi_{0}$, where $\phi_{0}$ satisfies

$P|\hat{\phi}_{n}(A,X)-\phi_{0}(A,X)|^{2}\xrightarrow{P}0,$ (2)

where $P$ means taking the expectation over $(X,A,Y)$ under the true model.

In general, most generalize linear models satisfy both assumptions under
regularity conditions. For example, the first assumption is automatically
satisfied for logistic regression. Suppose
$\pi(a|x,\theta_{0})=\frac{e^{a\theta_{0}^{T}x}}{1+e^{\theta_{0}^{T}x}}$, for
$a\in\\{0,1\\}$. The left hand side of the equation 1 becomes
$\sqrt{n}(\frac{e^{a\hat{\theta}_{n}^{T}x}}{1+e^{\hat{\theta}_{n}^{T}x}}-\frac{e^{a\theta_{0}^{T}x}}{1+e^{\theta_{0}^{T}x}})$.
By using the Taylor expansion for $\hat{\theta}_{n}$ at $\theta_{0}$, with
$a=1$ and $a=0$, the left hand side
$=\sqrt{n}(\hat{\theta}_{n}-\theta_{0})^{T}(2a-1)x\frac{e^{a\theta_{0}^{T}x}}{(1+e^{\theta_{0}^{T}x})^{2}}+o_{P}(1)$.
Thus
$\phi_{0}(a,x)=x(2a-1)\frac{e^{a\theta_{0}^{T}x}}{(1+e^{\theta_{0}^{T}x})^{2}}$
and
$\hat{\phi}_{n}(a,x)=x(2a-1)\frac{e^{a\hat{\theta}_{n}^{T}x}}{(1+e^{\hat{\theta}_{n}^{T}x})^{2}}$.
Assumption 1 is thus satisfied. And we can also show that the other part of
our assumptions hold by standard empirical process arguments. Note that in
randomized trial data, the propensity score will be fixed for binary treatment
$\pi(a|x)=0.5$ and $\hat{\theta}_{n}=0$.

###  2.2 Value Function

We evaluate the goodness of a decision rule $\hat{d}_{j,n}$ by estimating the
value function
$V(\hat{d}_{j,n})=E(Y|A=\hat{d}_{j,n}(X),X)=E(\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{P(A|X)})$
4. The value function can be interpreted as the inverse propensity score
weighted average outcome if the population were to follow the decision rule
$\hat{d}_{j,n}$. Specifically, we will calculate the value function based on
the test set data, represented by the subscript m:

$\hat{V}_{m}(\hat{d}_{j,n}(X_{i}))=\frac{\sum_{i=1}^{m}\frac{y_{i}1\\{A_{i}=\hat{d}_{j,n}(x_{i})\\}}{\hat{\pi}_{n}(A_{i}|X_{i})}}{\sum_{i=1}^{m}\frac{1\\{A_{i}=\hat{d}_{j,n}(x_{i})\\}}{\hat{\pi}_{n}(A_{i}|X_{i})}},$

where the propensity score $\hat{\pi}_{n}$ and the ITR $\hat{d}_{j,n}$ are
estimated from the training set and applied to the test set. The data
$(X_{i},A_{i},Y_{i})$, $i=1,\dots,m$ come from the test set. It is crucial to
ensure independence between training and test sets for the test set results to
effectively reflect the generalizability of the training set result.

###  2.3 Variance for the Value Function

######  Proposition 1.

$\begin{split}\sqrt{m}(\hat{V}_{m}(\hat{d}_{j,n})-V_{0}(\hat{d}_{j,n}))&=\mathbb{G}_{m}(\frac{(Y-V_{0}(\hat{d}_{j,n}))1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)})\\\
&-\sqrt{m/n}\sqrt{n}(\hat{\theta}_{n}-\theta_{0})^{T}E(\phi_{0}(A,X)\frac{(Y-V_{0}(\hat{d}_{j,n}))1\\{A=\hat{d}_{j,n}(x)\\}}{\pi^{2}_{0}(A|X)}),\\\
\end{split}$ (3)

where m denotes the empirical estimate based on the test set with $m$
individuals, $\mathbb{G}_{m}(U)=\sqrt{m}[\mathbb{P}_{m}(U)-P(U)]\rightarrow
N(0,Var(U))$ by standard empirical process arguments,
$\mathbb{P}_{m}(U)=m^{-1}\sum_{i=1}^{m}U_{i}$ is the empirical measure,
$P(U)=E(U)$ is the expectation taken over $U$, and
$\sqrt{n}(\hat{\theta}_{n}-\theta_{0})\rightarrow N(0,\Sigma_{0})$ by
Assumption 1.

Let
$U_{0}=\frac{(Y-V_{0}(\hat{d}_{j,n}))1\\{A=\hat{d}_{j,n}(x)\\}}{\pi_{0}(A|X)}$
and
$W_{0}=E(\phi_{0}(A,X)\frac{(Y-V_{0}(\hat{d}_{j,n}))1\\{A=\hat{d}_{j,n}(x)\\}}{\pi^{2}_{0}(A|X)})$.
Standard empirical process arguments yield
$\mathbb{G}_{m}(U_{0})=\sqrt{m}[\mathbb{P}_{n}(U_{0})-P(U_{0})]\rightarrow
N(0,Var(U_{0}))$. We can also show that the random variables
$\sqrt{m}(\hat{V}_{m}(\hat{d}_{j,n})-V_{0}(\hat{d}_{j,n}))$ converge in
distribution to a mean zero normal distribution:

$\sqrt{m}(\hat{V}_{m}(\hat{d}_{j,n})-V_{0}(\hat{d}_{j,n}))\rightarrow
N(0,Var(U_{0})+\frac{m}{n}W_{0}^{T}\Sigma_{0}W_{0}).$

In this context, as the value of $n$ approaches infinity, we assume the
quotient $m/n$ asymptotically approaches a finite limit, rather than diverging
towards infinity. Thus the expected value of the variance of a single ITR
$\hat{d}_{j,n}$ is:

$Var(\hat{V}_{m}(\hat{d}_{j,n}))=m^{-1}(Var(U_{0})+\frac{m}{n}W_{0}^{T}\Sigma_{0}W_{0})=m^{-2}\sum_{i=1}^{m}(\hat{U}_{i,j}-\bar{U}_{i,j})^{2}+\frac{1}{n}W_{0}^{T}n\hat{\Sigma}_{n}W_{0}.$

The variance $\sigma_{0}^{2}$ for the value function $V_{m}(d_{j})$ is then
estimated by:

$\hat{\sigma}^{2}_{m,j}=Var(\hat{V}_{m}(\hat{d}_{j}))=m^{-2}\sum_{i=1}^{m}(\hat{U}_{i,j}-\bar{U}_{i,j})^{2}+\frac{1}{n}\hat{W}_{j,m}^{T}n\hat{\Sigma}_{n}\hat{W}_{j,m}.$
(4)

Note that the two parts in the variance equation are dependent. When the
number of individuals in the test set is much less than that in the training
set, then we can ignore the second term
$\frac{1}{n}\hat{W}_{j,m}^{T}n\hat{\Sigma}_{n}\hat{W}_{j,m}$. When we use
randomized trial data, we don’t need to estimate the propensity score, and
thus $\Sigma_{n}=Var(\hat{\theta}_{n})=0$, and hence only the first term
remains. The estimated influence function for the training set and estimate
for $W_{0}$ are as follows:

$\displaystyle\hat{U}_{i,j}=(Y_{i}-\hat{V}_{n}(\hat{d}_{j,n}))\frac{1\\{A_{i}=\hat{d}_{j,n}\\}}{\hat{\pi}_{n}(A_{i}|X_{i})},\quad
i=1,...,n,$ (5)
$\displaystyle\hat{W}_{j,n,}=n^{-1}\sum_{i=1}^{n}[\frac{\hat{\phi}(A_{i},X_{i})\hat{U_{i,j}}}{\hat{\pi}_{n}(A_{i}|X_{i})}],$
(6)

where
$\hat{\phi}_{n}(A_{i},X_{i})=x_{i}(2a_{i}-1)\frac{e^{a_{i}\hat{\theta}_{n}^{T}x_{i}}}{(1+e^{\hat{\theta}_{n}^{T}x_{i}})^{2}}$,
$\bar{U}_{j,m}=m^{-1}\sum_{i=1}^{m}\hat{U}_{i,j}$ (j=1,2), and
$\hat{\Sigma}_{n}=Var(\hat{\theta}_{n})$ is estimated from the training set
based on the variance due to modeling the propensity score.

###  2.4 Comparison Between Value Functions

The distribution of the random variable for the difference between the two
value functions of two ITRs $\hat{d}_{1,n}$ and $\hat{d}_{2,n}$ can be shown
to converge to a normal distribution:

######  Proposition 2.

$\sqrt{m}(\hat{V}_{m}(\hat{d}_{1,n})-\hat{V}_{m}(\hat{d}_{2,n})-V_{0}(\hat{d}_{1,n})+V_{0}(\hat{d}_{2,n}))\rightarrow
N(0,T_{0}^{2}),$ (7)

where the empirical estimates of $T_{0}^{2}$ is:

$\hat{T}_{m}^{2}=m^{-1}\sum_{i=1}^{m}(\hat{U}_{i,1}-\hat{U}_{i,2}-\bar{U}_{m,1}+\bar{U}_{m,2})^{2}+m(\hat{W}_{m,1}-\hat{W}_{m,2})^{T}\hat{\Sigma}_{n}(\hat{W}_{m,1}-\hat{W}_{m,2}).$

Here, $\hat{\Sigma}_{n}$ is estimated from the training set based on both
variances due to multiple imputation and variance due to modeling the sample
data:

$\hat{\Sigma}_{n}=Var(\hat{\theta}_{n})=K^{-1}\sum_{k=1}^{K}VarCov(\hat{\theta}_{n,k})+(1+1/K)(K-1)^{-1}\sum_{k=1}^{K}(\hat{\theta}_{n,k}-\bar{\theta}_{n})^{2},$

where $\bar{\theta}_{n}=1/K\sum_{k=1}^{K}\hat{\theta}_{n,k}$.
$\hat{\Sigma}_{n}$ is the same across imputations. It is positive definite
because it is the sum of a positive definite matrix and a positive semi-
definite matrix.

Thus, the variance for the difference between the value function of two ITRs
is:

$Var(\hat{V}_{m}(\hat{d}_{1,n})-\hat{V}_{m}(\hat{d}_{2,n}))=m^{-2}\sum_{i=1}^{m}(\hat{U}_{i,1}-\hat{U}_{i,2}-\bar{U}_{m,1}+\bar{U}_{m,2})^{2}+(\hat{W}_{m,1}-\hat{W}_{m,2})^{T}\hat{\Sigma}_{n}(\hat{W}_{m,1}-\hat{W}_{m,2}).$
(8)

Under the null hypothesis, $V_{0}(\hat{d}_{1,n})-V_{0}(\hat{d}_{2,n})=0$. The
t-test statistics $t$ for this null hypothesis can be constructed as follows:

$t=\frac{\hat{V}_{m}(\hat{d}_{1,n})-\hat{V}_{m}(\hat{d}_{2,n})}{T_{m}/\sqrt{m}}\sim
N(0,1)$ (9)

The p-value for this t-test can be obtained by the test statistic $t$.

###  2.5 Multiple Imputation

In the presence of missing data, we extend our method to address the variance
of the value function in the case of multiple imputation. This extension is
based on the ideas from Chapter 2.3.2 from 17. Suppose we have K imputations.
Let $k=1,\dots,K$ be the index of the multiple-imputed data. Let
$\hat{V}_{m,k}(\hat{d}_{j})$ denote the value function for the $k^{th}$
imputed test set. Let $\hat{\sigma}_{m,j,k}=Var(\hat{V}_{m,k}(\hat{d}_{j}))$
be the estimator of the variance-covariance matrix for the estimated value
function for the $k^{th}$ imputed test set. Let
$\hat{T}^{2}_{m,k}=Var(\hat{V}_{m,k}(\hat{d}_{1})-\hat{V}_{m,k}(\hat{d}_{2}))$
denote the variance for the difference between two value functions for the
$k^{th}$ imputed test set. The subscript m indicates that the estimates are
obtained for the test set.

Let
$\tilde{V}_{m}(\hat{d}_{j})=K^{-1}\sum_{k=1}^{K}\hat{V}_{m,k}(\hat{d}_{j})$,
be value function estimate for decision rule $\hat{d}_{j}$. Let
$\sigma_{0,K}^{2}=Var_{K}(\hat{V}(\hat{d}_{j}))$ denote the variance
associated with the estimate for a single value function. The empirical
estimate of $\sigma_{0,K}^{2}$ is:

$\hat{\sigma}_{m,K}^{2}=(1+1/K)(K-1)^{-1}\sum_{k=1}^{K}(\hat{V}_{m,k}(\hat{d}_{j})-\tilde{V}_{m}(\hat{d}_{j}))^{2}+K^{-1}\sum_{k=1}^{K}\hat{\sigma}^{2}_{m,j,k}.$
(10)

For the paired test, Let
$T_{0,K}^{2}=Var_{K}(\hat{V}(\hat{d}_{1})-\hat{V}(\hat{d}_{2}))$ denote the
variance for the difference between value functions for decision rule
$\hat{d}_{1}$ and $\hat{d}_{2}$. The empirical estimate of $T_{0,K}^{2}$ is:

$\hat{T}_{m,K}^{2}=(1+1/K)(K-1)^{-1}\sum_{k=1}^{K}(\hat{V}_{m,k}(\hat{d}_{1})-\hat{V}_{m,k}(\hat{d}_{2})-\tilde{V}_{m}(\hat{d}_{1})+\tilde{V}_{m}(\hat{d}_{2}))^{2}+K^{-1}\sum_{k=1}^{K}\hat{T}^{2}_{m,k}.$
(11)

Under the null hypothesis, $V_{0,K}(\hat{d}_{1,n})-V_{0,K}(\hat{d}_{2,n})=0$.
The t-test statistics $t_{K}$ for this null hypothesis can be constructed as
follows:

$t_{K}=\frac{K^{-1}\sum_{k=1}^{K}(\hat{V}_{m,k}(\hat{d}_{1,n})-\hat{V}_{k,m}(\hat{d}_{2,n}))}{\sqrt{Var_{K}(\hat{V}(\hat{d}_{1})-\hat{V}(\hat{d}_{2}))}}\sim
N(0,1)$ (12)

The p-value for this t-test can be obtained from the test statistic $t_{K}$.

##  3 Simulation

Figure 1: The value function results of four simulation scenarios: (a) with
random treatment assignment and complete data, (b) with random treatment
assignment and data missing at random (c) with treatment assignment depending
on the covariates and complete data, (d) with treatment assignment depending
on the covariates and data missing at random. Within each scenario, we compare
the value functions of five different treatment regimes which are indicated in
the x-axis: observed treatment, treatment $A=0$ for all individuals, treatment
$A=1$ for all individuals, Q-learning optimal ITR, and D-learning optimal ITR.
For comparison purposes, we include the average outcome as a benchmark. The
true value function is presented with the grey line and the variation of the
true value function across the 10 replicates is indicated by the grey shaded
area.

Figure 2: The differences of value functions between selected pairs of the
following treatments: observed treatment, universal treatment $A=0$, universal
treatment $A=1$, Q-learning’s optimal ITR, and D-learning’s optimal ITR. The
results are presented under four simulation scenarios.

We evaluate and compare the value functions obtained from applying each of the
following treatment regimes to the population under two scenarios: (1) the
observed treatment, (2-3) the one-size-fits-all approach (one for each of the
binary treatments), (4) Q-learning ITR, (5) D-learning ITR. A convex penalty,
elastic net penalization 18, was used for each of the ITR models for variable
selection. We generate a total of 5000 observations for each simulated sample
with 70% for the training set and 30% for the test set. The number of
covariates included in the model is $p=20$. We perform 10 replicates for each
of the simulations.

The following data generation process is modified from the first simulation
scenario in 5. The covariates $X$ are uncorrelated and generated from a
multivariate normal distribution with mean zero and the covariance-variance
matrix is a diagonal of 1s. Let $A$ denote the binary treatment. In the model
implementation step, we will use $A=\\{0,1\\}$ for Q-learning and switch to
$A=\\{-1,1\\}$ for D-learning. Let $Y=\\{0,1\\}$ denote the binary outcome.
$\varepsilon\sim N(0,1)$ is the random error for generating the outcome $Y$.
Let $\beta_{0}=(\sqrt{6})^{-1}$, $\beta_{1}=\beta_{2}=0$,
$\beta_{j}=(2\sqrt{6})^{-1}$, $j=3,4,\dots,10$, $\beta_{11}=...=\beta_{p}=0$
be the coefficient for the main effect. Let
$(\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4},\gamma_{5},\dots,\gamma_{p})=(1,1,-1,1,-1,0,\dots,0)$
be the coefficients for treatment interactions. Let $Y\sim binomial$ with
probability equal to the expit of

$(\beta_{0}+\sum_{j=1}^{p}\beta_{j}X_{j})^{2}+(\gamma_{0}+\sum_{j=1}^{p}\gamma_{j}X_{j})A+\varepsilon.$

The generated data is designed to have small main effects
$(\beta_{0},\dots,\beta_{p})$, a large treatment effect $\gamma_{0}$, and
large treatment interaction effects $(\gamma_{1},\dots,\gamma_{p})$. To align
with the application example presented in section 4, we assume, without loss
of generality, that a smaller outcome Y and a smaller value function are
preferable. We consider four scenarios:

(a) Propensity score $\pi=0.5$ (randomized trial data)

(b) Propensity score $\pi=0.5$ (randomized trial data with data missing at
random)

(c) Propensity score $logit(\pi)=0.75X_{1}-0.75X_{2}$ (observational data)

(d) Propensity score $logit(\pi)=0.75X_{1}-0.75X_{2}$ (observational data with
data missing at random)

In scenarios where data are missing at random, the missingness mechanism is
implemented by assigning a conditional probability of missingness to the
variables $X_{2},\ldots,X_{p}$. This probability is contingent upon the value
of $X_{1}$: specifically, there is a 15% probability of missingness when
$X_{1}>0$, and a notably lower probability of 10% in cases where this
condition is not met. Under each scenario, the true value function is
estimated based on a simulated large test set with 10000 individuals. The data
generation process for this test set is the same as for the training set under
the same scenario. The true variance for the value function is
$\hat{\tau}_{j}^{2}=\frac{1}{L}\sum_{l=1}^{L}(V_{l}(\hat{d}_{j})-V_{*}(\hat{d}_{j}))^{2}$.
$V_{l}(\hat{d}_{j})$ is the estimated value function based on the test set $l$
where we plug in the estimated decision rule $\hat{d}_{j}$ and the estimated
propensity score model from the training set. $V_{*}(\hat{d}_{j})$ is the
estimated value function based on the giant test set, and is based on the true
propensity score parameters and the decision rule $\hat{d}_{j}$ estimated from
the training set. The true variance for the difference between value functions
is constructed similarly.

Figure 1 shows the value functions across these four scenarios. We observe
that the value function for ITRs from Q-learning and D-learning are the lowest
among all methods. The value function estimates associated with the optimal
ITRs estimated using Q-learning and D-learning are close to the true value
function. The estimated variance associated with these two value functions is
also close to the true variance around the true value function. Since
Q-learning shows the best value function and lowest optimal treatment
misclassification rate overall (Table 1), we focus on comparing Q-learning
with the others in the pairwise comparison as shown in Figure 2. We observe
that Q-learning is significantly better than the one-size-fits-all treatments
($A=0$ and $A=1$) in all four scenarios. This is expected because we simulate
the data with a large prescriptive treatment effect. We observe that the
difference is large between value functions for assigning treatment $A=0$ and
$A=1$ to all individuals. This is expected because we simulate the data with a
large treatment effect. For observational data, when treatment assignments are
nonrandom and probabilities need to be estimated by propensity scores, value
function estimates have larger variances compared to the value function
estimates for randomized trial data (Table 1). For datasets with missing
values, the performance remains comparable to that of complete datasets when
the missing data is random and is addressed through multiple imputation using
missForest. The results shown here are the test set results. The training set
results are similar to the test set results but with a smaller variance. The
simulation results illustrate the effectiveness of our approach in discerning
disparities between two ITRs while maintaining reasonable variance estimates
in both observational data and randomized trial data, with or without missing
values. Additionally, it highlights a constraint of our method that the
accuracy of the estimate’s variance relies on the precision of propensity
score estimation.

Table 1: Estimated value functions and treatment misclassification rate (MC) of different treatment regimes and under four scenarios for simulation data. | Scenario 1 | Scenario 2 | Scenario 3 | Scenario 4
---|---|---|---|---
| Mean(SD) | MC | Mean(SD) | MC | Mean(SD) | MC | Mean(SD) | MC
Obs | 0.631(0.026) | 0.509 | 0.631(0.026) | 0.509 | 0.631(0.028) | 0.581 | 0.631(0.027) | 0.581
All0 | 0.567(0.018) | 0.235 | 0.567(0.018) | 0.235 | 0.561(0.020) | 0.235 | 0.562(0.019) | 0.235
All1 | 0.696(0.017) | 0.765 | 0.696(0.017) | 0.765 | 0.699(0.020) | 0.765 | 0.700(0.019) | 0.765
Q-Learning | 0.534(0.019) | 0.064 | 0.535(0.019) | 0.065 | 0.529(0.022) | 0.063 | 0.532(0.021) | 0.066
D-Learning | 0.536(0.019) | 0.087 | 0.537(0.019) | 0.089 | 0.534(0.022) | 0.094 | 0.535(0.021) | 0.095

##  4 Application

###  4.1 CHNS: Data Analysis

Hypertension, also known as high blood pressure, is a chronic medical
condition that affects over one billion adults worldwide 19. It is a major
risk factor for cardiovascular disease, which is the leading cause of death
globally 20. People with hypertension are more likely to develop kidney
disease 21, vision problems 22, and cognitive impairment 23. Hypertension can
also affect quality of life and increase annual medical costs 24. Over the
past two to three decades, China has experienced a significant increase in the
prevalence of hypertension (from 20.8% in 2004 to 29.6% in 2010 and 24.7% in
2018) due to increased life expectancy and lifestyle changes 25, 26. As
hypertension has become a significant burden on China’s population health, it
is important to study this condition and develop effective prevention and
treatment strategies. While medication proves effective in managing
hypertension, the general adherence to anti-hypertensive medication remains
low. This can be attributed to factors such as a low reimbursement ratio, high
daily medical costs, and access to healthcare 27. In addition to taking
medications, adopting non-pharmacologic interventions such as engaging in
regular physical activity, reducing salt intake, and following a balanced diet
rich in fruits and vegetables can help to reduce high blood pressure 28.
However, there is potential heterogeneity in which individuals may see the
most benefit from behavioral interventions.

Our overall objective is to identify the most effective non-pharmacologic
treatments for subpopulations of individuals, considering the challenges of
adhering to a comprehensive healthy lifestyle based on public health
recommendations. As a first step, we examine the potential variability in the
effectiveness of exceeding 5 hours of weekly moderate-to-vigorous physical
activity (MVPA) in reducing the risk of hypertension in the CHNS population.
We seek to distinguish between participants who benefit from engaging in MVPA
for more than 5 hours per week and those who do not. For participants who see
no hypertension benefit from increasing their physical activity, it is
possible to identify other potential interventions to consider, such as diet
modification. Understanding which interventions confer optimal outcomes
empowers individuals to identify the most effective and efficient behavioral
treatments and their associated benefits. This knowledge can be highly
beneficial, as individuals are more motivated to engage in behaviors when they
are aware of their positive outcomes 29. To address this research question, we
tested whether a personalized physical intervention strategy would improve
population-level health outcomes compared to a one-size-fits-all intervention.

We used data from the China Health and Nutrition Survey (CHNS) 30, a
population-based, observational data set consisting of high-quality data on
diet (3-day repeated 24-hour recalls) and physical activity (detailed 7-day
recall) collected from individuals in China using detailed recall instruments.
Our analysis focused on the study year 2009 with 8320 adults. Our outcome,
hypertension, was defined 111Anti-hypertensive medication is included in the
model in our analysis, so blood pressure measurements here are not adjusted
for medication. Participants who were taking medication may exhibit varying
responses to changes in physical activity. Consequently, we have included
medication as a covariate that has the potential to influence how individuals
respond to engaging in more than 5 hours of physical activity per week. as a
systolic blood pressure $\geq$ 130 mmHg or diastolic blood pressure $\geq$ 80
mmHg. Individuals with missing systolic or diastolic blood pressure
measurements or who were missing more than 10% of their covariates were
excluded from our study, resulting in a final analytic sample of 5241
individuals, 60% of whom had hypertension. In the CHNS dataset, physical
activity encompassing leisure, occupational, transportation, and domestic
activities, were collected via a comprehensive survey recall instrument 31.
Our study focused on weekly moderate-to-vigorous physical activity (MVPA),
defined as activities with a minimum of 3 METS (Metabolic equivalent of task)
32 222The METS, or Metabolic equivalent of task, is a unit of measurement for
physical activity. One MET is equivalent to a person’s oxygen consumption at a
rate of 3.5 milliliters per kilogram per minute. 33. Our sample population has
a high MVPA level, with a median of 15 hours per week 31. Targeting the 5-hour
MVPA weekly guideline 2, 1, 31, we dichotomized the treatment into a binary
variable $A$: over 5 hours ($A=1$) and 5 or fewer hours ($A=0$) of MVPA per
week. In the CHNS data, 66% of the sample population engaged in more than 5
hours of MVPA weekly. With insights from biological understanding and
published literature 31, we identified twenty risk factors for hypertension.
These factors also have the potential to influence the relationship between
physical activity and hypertension, including: age, gender, caloric intake,
BMI, education, smoking status, sodium intake, potassium intake, alcohol
consumption, hypertension medication, household income, province, and
urbanization index. Less than 15% of our sample had missing covariates,
handled using missForest package in R for multiple imputation 34, 35. Table 2
provides further details about these factors based on one imputed dataset. The
other imputed datasets exhibit similar distributions.

Table 2: Descriptive Table for CHNS Data.

| Training Data | Test Data | Overall
---|---|---|---
| (N=3667) | (N=1574) | (N=5241)
Age, mean (SD) | 48.7 (10.4) | 48.8 (10.2) | 48.7 (10.3)
Gender, n (%) |  |  |
Male | 1695 (46.2%) | 746 (47.4%) | 2441 (46.6%)
Female | 1972 (53.8%) | 828 (52.6%) | 2800 (53.4%)
Calorie Intake (kcal/kg), mean (SD) | 36.7 (11.8) | 36.7 (11.8) | 36.7 (11.8)
BMI, mean (SD) | 23.3 (3.18) | 23.4 (3.31) | 23.3 (3.22)
Education, mean (SD) | 1.56 (0.710) | 1.54 (0.723) | 1.56 (0.714)
Log-transformed Household Income, mean (SD) | 2.51 (0.853) | 2.50 (0.849) | 2.51 (0.852)
Current Smoking, n (%) | 1234 (33.7%) | 510 (32.4%) | 1744 (33.3%)
Sodium (mg/day), mean (SD) | 4560 (2270) | 4600 (2250) | 4570 (2260)
Potassium (mg/day), mean (SD) | 1750 (673) | 1720 (657) | 1740 (668)
Alcohol, n (%) | 1262 (34.4%) | 572 (36.3%) | 1834 (35.0%)
MVPA (hours/week) $\dagger$, mean (SD) | 34.2 (42.0) | 36.3 (43.3) | 34.8 (42.4)
Anti-Hypertensive Medication, n (%) | 291 (7.9%) | 144 (9.1%) | 435 (8.3%)
Hypertension, n (%) | 2184 (59.6%) | 940 (59.7%) | 3124 (59.6%)
Province, n (%) |  |  |
21 – Liaoning | 400 (10.9%) | 170 (10.8%) | 570 (10.9%)
23 – Heilongjiang | 464 (12.7%) | 199 (12.6%) | 663 (12.7%)
32 – Jiangsu | 424 (11.6%) | 159 (10.1%) | 583 (11.1%)
37 – Shandong | 422 (11.5%) | 200 (12.7%) | 622 (11.9%)
41 – Henan | 335 (9.1%) | 148 (9.4%) | 483 (9.2%)
42 – Hubei | 367 (10.0%) | 172 (10.9%) | 539 (10.3%)
43 – Hunan | 444 (12.1%) | 183 (11.6%) | 627 (12.0%)
45 – Guangxi | 415 (11.3%) | 156 (9.9%) | 571 (10.9%)
52 – Guizhou | 396 (10.8%) | 187 (11.9%) | 583 (11.1%)
Urbanization Index, mean (SD) | 66.1 (19.0) | 65.3 (18.7) | 65.9 (18.9)

  * •

Source: China Health and Nutrition Survey (CHNS).

  * $\dagger$

MVPA includes occupational, domestic, transportation and leisure activity.

Figure 3: The value function results of six different treatment regimes:
observed treatment, treatment $A=0$ for all individuals, treatment $A=1$ for
all individuals, Q-learning optimal ITR, and D-learning optimal ITR. For
comparison purposes, we include the average outcome as a benchmark. (a) shows
the results on original data, and (b)(c)(d) shows the three simulated data
with three levels of treatment effect modification $\delta=1$, $\delta=2$,
$\delta=10$, respectively.

Figure 4: The value function differences among selected pairs of treatment
regimes. Specifically, the Q-learning optimal ITR comparison with the observed
treatment, treatment $A=0$ for all individuals, treatment $A=1$ for all
individuals, and D-learning optimal ITR are shown. (a) shows the results on
original data, and (b)(c)(d) shows the three simulated data with three levels
of treatment effect modification $\delta=1$, $\delta=2$, $\delta=10$,
respectively.

The data was split into a 70% training set and a 30% test set, with missForest
imputation performed 10 times on each set to avoid dependency. We obtained
Q-learning optimal ITRs 4 and D-learnings optimal ITR 5 for each of the
training sets using a logistic regression model with elastic-net penalization.
We then compared the value function of each treatment rule on the test set,
including the observed treatment, assigning all individuals to treatment $A=1$
or $A=0$, optimal ITRs derived from Q-learning and D-learning. The value
functions of these treatment rules were estimated using the common estimator
$V$ and the differences between value functions were directly obtained by
subtracting the value functions and then averaged over 10 imputed data sets.
The variance of the value function and the variance of the difference between
the value functions were estimated by our method, which accounts for both
modeling and multiple imputation variance.

Table 3: Estimated value functions of different treatment regimes for CHNS data. | Training Set | Test Set
---|---|---
| Mean | SD | Mean | SD
Obs | 0.602 | 0.021 | $0.595$ | 0.031
AllLowPA | 0.603 | 0.018 | $0.591$ | 0.026
AllHighPA | 0.600 | 0.011 | $0.598$ | 0.017
Q-Learning | 0.568 | 0.013 | $0.574$ | 0.017
D-Learning | 0.576 | 0.016 | $0.598$ | 0.023

Table 4: Estimated value differences and p-values of selected comparisons between treatment regimes for CHNS data. | Training Set | Test Set
---|---|---
Difference | Mean | SD | P-Value | Mean | SD | P-Value
AllLowPA vs. Obs | 0.001 | 0.011 | 0.924 | -0.004 | 0.017 | 0.827
AllHighPA vs. Obs | -0.001 | 0.018 | 0.953 | 0.003 | 0.026 | 0.897
AllHighPA vs. AllLowPA | -0.002 | 0.021 | 0.921 | 0.007 | 0.032 | 0.821
Q vs. Obs | -0.033 | 0.017 | 0.046 | -0.021 | 0.026 | 0.429
Q vs. AllLowPA | -0.034 | 0.017 | 0.040 | -0.017 | 0.026 | 0.519
Q vs. AllHighPA | -0.032 | 0.013 | 0.014 | -0.024 | 0.018 | 0.176
Q vs. D | -0.008 | 0.013 | 0.562 | -0.023 | 0.020 | 0.239

Our goal is to compare personalized and uniform interventions in reducing
hypertension risk by comparing the value functions in test sets for (1)
observed physical activity, under the assumption that individuals continue
their existing routines; (2) assigning all individuals to treatment $A=1$,
assuming all individuals are above the 5 hours weekly MVPA recommendation; (3)
assigning all individuals to treatment $A=0$, assuming all individuals are at
or below the 5 hours weekly MVPA recommendation; and (4) optimal ITRs derived
from various methods. Figure 3(a) shows the value function and variances for
those treatment regimes. The Q-learning value function is the lowest. This
means that if the population follows the Q-learning ITR, the expected risk of
hypertension in the population is lower than that if the population follows
the other treatment assignment, although the differences between the value
functions are small. Figure 4(a) shows the pairwise difference of the value
functions between treatment regimes that we are interested in. Our focus is on
determining the necessity of recommending over 5 hours of weekly MVPA to
everyone or specific subgroups. Thus we compare the best ITR, Q-learning ITR,
defined by the lowest value function estimate in the training set (Table 3)
with the one-size-fits-all treatment $A=1$ for the test set. Based on the
t-test result, the value functions are not significantly different, with a
p-value of 0.18 for the Q-learning personalized treatment rule compared to
assigning all individuals to treatment $A=1$ (Table 4). This suggests that a
personalized approach does not significantly outperform the population-level
physical activity recommendations for doing physical activity for more than
five hours per week. Although the training data results show that the
Q-learning has a significantly lower value function than the best one-size-
fits-all treatment $A=1$, the results shown here are test set results. The
test set results hold particular significance as they underscore the model’s
generalizability. They offer insight into the performance expected of the
trained model when confronted with new and unseen data.

Our findings suggest that adjustments in physical activity for a subset of
individuals alone do not yield a substantial reduction in hypertension risk at
the population level. Alternative non-pharmacologic interventions may
demonstrate greater efficacy in mitigating hypertension risk across specific
subpopulations. This variance in response may be influenced by the interaction
between multiple factors, such as sodium intake and physical activity. For
instance, individuals with high sodium consumption may not benefit
significantly from increased physical activity alone. In contrast, sodium
intake reduction could potentially offer a more effective strategy for
lowering hypertension risk. Another possible explanation is that much physical
activity in the CHNS data comprised working, as opposed to leisure,
activities. Published data suggest that occupational physical activity may not
be as protective against hypertension 36. The complete pairwise comparison of
all the comparisons is depicted in the appendix Figure 6.

###  4.2 CHNS: Augmented Analysis

These differences across all the value functions in our original CHNS analysis
could truly be this small in magnitude. However, it is important to note that
these models only consider a small subset of factors that could impact how
changes in physical activity can impact hypertension. Hypertension is complex,
resulting from an interplay between individual behavior, environment, and
additional factors like metabolism and genetics. These factors might help
provide additional important insight and better identify underlying
heterogeneity in the impact of physical activity on hypertension not captured
in our original analysis. Therefore, we wanted to investigate our method’s
performance in the presence of prescriptive treatment effect interaction,
which corresponds to these additional factors, or combinations of these
additional factors with our existing factors, that explain more of the
variation in how an individual’s hypertension risk would be impacted by
changes in physical activity. Therefore, an additional analysis was carried
out on simulated datasets with modified prescriptive treatment effect
interaction terms.

The outcome variables in these data sets were generated by the fitted
Q-learning outcome regression model with the treatment effect interaction
terms (excluding the treatment main effect term) multiplied by a factor of
$\delta\geq 1$. We investigate three options of $\delta$ for this analysis:
$\delta=1$, $\delta=2$, and $\delta=10$, representing the estimated original,
double, and tenfold treatment effect, respectively. In each analysis,
hypertension outcomes are replaced by generated outcomes from the modified
Q-learning outcome regression model. This augmented analysis examines our
method’s performance when the outcome generative process is known and how the
performance varies as the true prescriptive treatment effect enlarges.

As $\delta$ increases, from $\delta=1$, to $\delta=2$, to $\delta=10$ (Figure
3 (b)-(d)), the value function differences between the Q-learning strategy and
other treatment rules increase. Figure 4 shows the difference in magnitude
between the pairwise differences in value function results between the
Q-learning treatment rules and both (1) assigning all individuals to either
treatment $A=1$ and (2) treatment $A=0$ across our original results. The
complete pairwise comparison of all the comparisons is depicted in the
appendix figure 6. In all these instances, the Q-learning models deliver the
best ITRs, with values notably lower, suggesting a lower risk of hypertension
in the treatment-effect-augmented population especially when compared to
population interventions when all individuals are assigned to low physical
activity or all individuals are assigned to high physical activity. When
$\delta=1$, we observe that the pattern of value functions bears resemblance
to that of the original data across various methods. As the prescriptive
treatment effect becomes increasingly evident, both Q-learning and D-learning
demonstrate progressively superior performance compared to the ’one-size-fits-
all’ approaches. This underscores the importance of personalizing over
population-level interventions when certain factors significantly affect
individual responses to treatment. Applying these methods to existing rich,
detailed observational data can provide an important insight into factors that
might distinguish these subgroups of responders. Future work will build on
these existing models to better incorporate additional high-dimensional data.
However, these results suggest that if a large heterogeneous effect truly
exists, our method can successfully capture it.

##  5 Discussion

In this paper, we propose an innovative t-test based approach that can
directly compare the value functions of any two treatment regimes. The
validity of the approach follows from the asymptotic normality of the standard
value function for the estimated treatment regime and the asymptotic normality
of the propensity score model parameters. Our method provides valid estimates
for (1) the variance for a value function for a single treatment regime, (2)
the variance for the difference between two value functions of two treatment
regimes, (3) the p-value of the t-test for the significance of the difference
between two value functions, and (4) the application of these estimations in
scenarios involving multiple imputations. This method maintains simplicity in
variance calculation and is computationally more efficient than the bootstrap
method, especially when multiple imputation is used for handling missing data.
Through simulation studies and the data application example, we demonstrate
the performance and the ease of implementation of our method in different
scenarios. Additionally, this method enables the evaluation and comparison of
ITR effectiveness using abundant observational cohort data when clinical trial
data are not available. This comparison is crucial as individuals seek
behavioral modification guidance amidst numerous recommendations, and
emphasizing targeting strategies could increase the likelihood of effective
changes.

Nevertheless, it is important to acknowledge a limitation inherent in our
methodology – its reliance on the presupposition that the propensity score
model holds true. In observational data, extreme propensity scores can
influence the estimates of the value function as well as the variance
estimates of the value function. For future work, it will be interesting to
use more robust estimators to alleviate the impact of potential
misspecification of the propensity score model. Another extension of this
method is to multiple-stage ITR comparison. Expanding our method to
simultaneously compare more than two ITRs would also present an interesting
avenue for further exploration.

To conclude, our method offers a convenient approach to comparing two
treatment regimes directly. Furthermore, it is suitable for observational data
and randomized trial data and has the ability to incorporate multiple
imputation for missing data.

## Acknowledgments

This work was supported by the NIH, Eunice Kennedy Shriver National Institute
of Child Health and Human Development (R01 HD30880), and the National
Institute on Aging (R01AG065357). This research uses data from China Health
and Nutrition Survey (CHNS). We are grateful to research grant funding from
the National Institute for Health (NIH), the Eunice Kennedy Shriver National
Institute of Child Health and Human Development (NICHD) for R01 HD30880 and
R01 HD38700, National Institute on Aging (NIA) for R01 AG065357, National
Institute of Diabetes and Digestive and Kidney Diseases (NIDDK) for R01
DK104371 and P30 DK056350, National Heart, Lung, and Blood Institute (NHLBI)
for R01 HL108427, the NIH Fogarty grant D43 TW009077, the Carolina Population
Center for P2C HD050924 and P30 AG066615 since 1989, and the China-Japan
Friendship Hospital, Ministry of Health for support for CHNS 2009, Chinese
National Human Genome Center at Shanghai since 2009, and Beijing Municipal
Center for Disease Prevention and Control since 2011. We thank the National
Institute for Nutrition and Health, China Center for Disease Control and
Prevention, Beijing Municipal Center for Disease Control and Prevention, and
the Chinese National Human Genome Center at Shanghai. In addition, Minxin Lu
was supported by NC TraCS collaboration on "Optimizing weight status based on
potentially modifiable risk factors". Both Minxin Lu and Michael Kosorok were
funded in part by grant UM1 TR004406 from the National Center for Advancing
Translational Sciences. We thank Matthew Christopher Brown and Lina Maria
Montoya for their help and support for the project. The content is solely the
responsibility of the authors and does not necessarily represent the official
views of the NIH.

### Author contributions

Minxin Lu, Annie Green Howard, Penny Gordon-Larsen, Katie A. Meyer, Shufa Du,
and Michael R. Kosorok contributed to the study conception and design. Data
preparation was performed by Annie Green Howard and Hsiao-Chuan Tien. Data
collection was performed by Huijin Wang and Bing Zhang. Statistical analysis
was performed by Minxin Lu and Michael R. Kosorok. The first draft of the
manuscript was written by Minxin Lu and all authors commented on previous
versions of the manuscript. All authors read and approved the final
manuscript.

### Financial disclosure

None reported.

### Conflict of interest

The authors declare no potential conflict of interests.

## References

  * Piercy et al. 2018 Katrina L Piercy, Richard P Troiano, Rachel M Ballard, Susan A Carlson, Janet E Fulton, Deborah A Galuska, Stephanie M George, and Richard D Olson.  The physical activity guidelines for americans.  _Jama_ , 320(19):2020–2028, 2018.
  * Olson et al. 2019 Richard D Olson, Katrina L Piercy, Richard P Troiano, Rachel M Ballard, Janet E Fulton, Deborah A Galuska, Shellie Y Pfohl, Alison Vaux-Bjerke, Julia B Quam, Stephanie M George, Kyle Sprow, Susan A Carlson, Eric T Hyde, and Kate Olscamp.  _Physical Activity Guidelines for Americans 2nd edition_.  The U.S. Department of Health and Human Services, https://health.gov/sites/default/files/2019-09/Physical_Activity_Guidelines_2nd_edition.pdf, 2019.
  * Bouchard and Rankinen 2001 Claude Bouchard and Tuomo Rankinen.  Individual differences in response to regular physical activity.  _Medicine & Science in Sports & Exercise_, 33(6):S446–S451, 2001.
  * Qian and Murphy 2011 Min Qian and Susan A Murphy.  Performance guarantees for individualized treatment rules.  _Annals of statistics_ , 39(2):1180, 2011.
  * Tian et al. 2014 Lu Tian, Ash A Alizadeh, Andrew J Gentles, and Robert Tibshirani.  A simple method for estimating interactions between a treatment and a large number of covariates.  _Journal of the American Statistical Association_ , 109(508):1517–1532, 2014.
  * Sterne et al. 2009 Jonathan AC Sterne, Ian R White, John B Carlin, Michael Spratt, Patrick Royston, Michael G Kenward, Angela M Wood, and James R Carpenter.  Multiple imputation for missing data in epidemiological and clinical research: potential and pitfalls.  _Bmj_ , 338, 2009.
  * Jiang et al. 2021 Xiaotong Jiang, Amanda E Nelson, Rebecca J Cleveland, Daniel P Beavers, Todd A Schwartz, Liubov Arbeeva, Carolina Alvarez, Leigh F Callahan, Stephen Messier, Richard Loeser, et al.  Precision medicine approach to develop and internally validate optimal exercise and weight-loss treatments for overweight and obese adults with knee osteoarthritis: data from a single-center randomized trial.  _Arthritis care & research_, 73(5):693–701, 2021.
  * Cui et al. 2017 Yifan Cui, Ruoqing Zhu, and Michael Kosorok.  Tree based weighted learning for estimating individualized treatment rules with censored data.  _Electronic journal of statistics_ , 11(2):3927, 2017.
  * Zhao et al. 2012 Yingqi Zhao, Donglin Zeng, A John Rush, and Michael R Kosorok.  Estimating individualized treatment rules using outcome weighted learning.  _Journal of the American Statistical Association_ , 107(499):1106–1118, 2012.
  * Efron 1992 Bradley Efron.  Bootstrap methods: another look at the jackknife.  In _Breakthroughs in statistics: Methodology and distribution_ , pages 569–593. Springer, 1992.
  * Shi et al. 2022 Chengchun Shi, Sheng Zhang, Wenbin Lu, and Rui Song.  Statistical inference of the value function for reinforcement learning in infinite-horizon settings.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 84(3):765–793, 2022.
  * Bates et al. 2021 Stephen Bates, Trevor Hastie, and Robert Tibshirani.  Cross-validation: what does it estimate and how well does it do it?  _arXiv preprint arXiv:2104.00673_ , 2021.
  * Chakraborty et al. 2010 Bibhas Chakraborty, Susan Murphy, and Victor Strecher.  Inference for non-regular parameters in optimal dynamic treatment regimes.  _Statistical methods in medical research_ , 19(3):317–343, 2010.
  * Laber and Murphy 2011 Eric B Laber and Susan A Murphy.  Adaptive confidence intervals for the test error in classification.  _Journal of the American Statistical Association_ , 106(495):904–913, 2011.
  * Chakraborty et al. 2013 Bibhas Chakraborty, Eric B Laber, and Yingqi Zhao.  Inference for optimal dynamic treatment regimes using an adaptive m-out-of-n bootstrap scheme.  _Biometrics_ , 69(3):714–723, 2013.
  * Austin 2011 Peter C Austin.  An introduction to propensity score methods for reducing the effects of confounding in observational studies.  _Multivariate behavioral research_ , 46(3):399–424, 2011.
  * Van Buuren 2018 Stef Van Buuren.  _Flexible imputation of missing data_.  CRC press, 2018.
  * Zou and Hastie 2005 Hui Zou and Trevor Hastie.  Regularization and variable selection via the elastic net.  _Journal of the royal statistical society: series B (statistical methodology)_ , 67(2):301–320, 2005.
  * Bloch 2016 Michael J Bloch.  Worldwide prevalence of hypertension exceeds 1.3 billion.  _Journal of the American Society of Hypertension: JASH_ , 10(10):753–754, 2016.
  * Bromfield and Muntner 2013 Samantha Bromfield and Paul Muntner.  High blood pressure: the leading global burden of disease risk factor and the need for worldwide prevention programs.  _Current hypertension reports_ , 15:134–136, 2013.
  * Weldegiorgis and Woodward 2020 Misghina Weldegiorgis and Mark Woodward.  The impact of hypertension on chronic kidney disease and end-stage renal disease is greater in men than women: a systematic review and meta-analysis.  _BMC nephrology_ , 21(1):1–9, 2020.
  * Bhargava et al. 2012 M Bhargava, MK Ikram, and Tien Yin Wong.  How does hypertension affect your eyes?  _Journal of human hypertension_ , 26(2):71–83, 2012.
  * Kilander et al. 1998 Lena Kilander, Hakan Nyman, Merike Boberg, Lennart Hansson, and Hans Lithell.  Hypertension is related to cognitive impairment: a 20-year follow-up of 999 men.  _Hypertension_ , 31(3):780–786, 1998.
  * Wang et al. 2017 Guijing Wang, Xilin Zhou, Xiaohui Zhuo, and Ping Zhang.  Annual total medical expenditures associated with hypertension by diabetes status in us adults.  _American journal of preventive medicine_ , 53(6):S182–S189, 2017.
  * Zhang et al. 2023 Mei Zhang, Yu Shi, Bin Zhou, Zhengjing Huang, Zhenping Zhao, Chun Li, Xiao Zhang, Guiyuan Han, Ke Peng, Xinhua Li, et al.  Prevalence, awareness, treatment, and control of hypertension in china, 2004-18: findings from six rounds of a national survey.  _bmj_ , 380, 2023.
  * Wang et al. 2023 Ji-Guang Wang, Wei Zhang, Yan Li, and Lisheng Liu.  Hypertension in china: epidemiology and treatment initiatives.  _Nature Reviews Cardiology_ , pages 1–15, 2023.
  * Cui et al. 2020 Bin Cui, Zhaohui Dong, Mengmeng Zhao, Shanshan Li, Hua Xiao, Zhitao Liu, and Xiaowei Yan.  Analysis of adherence to antihypertensive drugs in chinese patients with hypertension: a retrospective analysis using the china health insurance association database.  _Patient preference and adherence_ , pages 1195–1204, 2020.
  * Appel 2003 Lawrence J Appel.  Lifestyle modification as a means to prevent and treat high blood pressure.  _Journal of the American Society of Nephrology_ , 14(suppl 2):S99–S102, 2003.
  * Cane et al. 2012 James Cane, Denise O’Connor, and Susan Michie.  Validation of the theoretical domains framework for use in behaviour change and implementation research.  _Implementation science_ , 7:1–17, 2012.
  * Zhang et al. 2014 Bing Zhang, FY Zhai, SF Du, and Barry M Popkin.  The c hina h ealth and n utrition s urvey, 1989–2011.  _Obesity reviews_ , 15:2–7, 2014.
  * Ng et al. 2014 Shu Wen Ng, A-G Howard, HJ Wang, Chang Su, and Bing Zhang.  The physical activity transition among adults in c hina: 1991–2011.  _Obesity Reviews_ , 15:27–36, 2014.
  * Ainsworth et al. 2011 Barbara E Ainsworth, William L Haskell, Stephen D Herrmann, Nathanael Meckes, David R Bassett Jr, Catrine Tudor-Locke, Jennifer L Greer, Jesse Vezina, Melicia C Whitt-Glover, and Arthur S Leon.  2011 compendium of physical activities: a second update of codes and met values.  _Medicine & science in sports & exercise_, 43(8):1575–1581, 2011.
  * Sylvia et al. 2014 Louisa G Sylvia, Emily E Bernstein, Jane L Hubbard, Leigh Keating, and Ellen J Anderson.  A practical guide to measuring physical activity.  _Journal of the Academy of Nutrition and Dietetics_ , 114(2):199, 2014.
  * Stekhoven 2022 Daniel J. Stekhoven.  _missForest: Nonparametric Missing Value Imputation using Random Forest_.  CRAN, https://cran.r-project.org/package=missForest, 2022.  R package version 1.5.
  * Stekhoven and Buehlmann 2012 Daniel J. Stekhoven and Peter Buehlmann.  Missforest - non-parametric missing value imputation for mixed-type data.  _Bioinformatics_ , 28(1):112–118, 2012.
  * Huai et al. 2013 Pengcheng Huai, Huanmiao Xun, Kathleen Heather Reilly, Yiguan Wang, Wei Ma, and Bo Xi.  Physical activity and risk of hypertension: a meta-analysis of prospective cohort studies.  _Hypertension_ , 62(6):1021–1026, 2013.

##  Appendix A Proofs

Proof of Proposition 1

###### Proof.

$\sqrt{m}(\hat{V}_{m}(\hat{d}_{j,n})-V_{0}(\hat{d}_{j,n}))=$

$\sqrt{m}(\frac{\sum_{i=1}^{m}\frac{y_{i}1\\{a_{i}=\hat{d}_{j,n}(x_{i})\\}}{\hat{\pi}_{n}(a_{i}|x_{i})}-mE(\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)})}{\sum_{i=1}^{m}\frac{1\\{a_{i}=\hat{d}_{j,n}(x_{i})\\}}{\hat{\pi}_{n}(a_{i}|x_{i})}})$

$-\sqrt{m}(\frac{E[\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X))}][\sum_{i=1}^{n}\frac{1\\{a_{i}=\hat{d}_{j,n}(x_{i})\\}}{\hat{\pi}_{n}(a_{i}|x_{i})}-mE(\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)})]}{E[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X))}][\sum_{i=1}^{n}\frac{1\\{a_{i}=\hat{d}_{j,n}(x_{i})\\}}{\hat{\pi}_{n}(a_{i}|x_{i})}]})$

$+\sqrt{m}\left(\frac{E[\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X))}]-E[\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X))}]}{E[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X))}]}\right)$

$-\sqrt{m}(\frac{E[\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X))}][E[\frac{1\\{A=\hat{d}_{j,n}(X\\}}{\hat{\pi}_{n}(A|X)}-E(\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)})]}{E[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X))}]E[\frac{1\\{A_{=}\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)}]})$

$=B_{1,m}-B_{2,m}+B_{3,m}-B_{4,m},$

where

$B_{1,m}=\sqrt{m}(\frac{\sum_{i=1}^{m}\frac{y_{i}1\\{a_{i}=\hat{d}_{j,n}(x_{i})\\}}{\hat{\pi}_{n}(a_{i}|x_{i})}-mE(\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)})}{\sum_{i=1}^{m}\frac{1\\{a_{i}=\hat{d}_{j,n}(x_{i})\\}}{\hat{\pi}_{n}(a_{i}|x_{i})}})$

$=\frac{\mathbb{G}_{n}(\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)})}{\mathbb{P}_{n}[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)}]}=\mathbb{G}_{n}[\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)}]+o_{P}(1),$

where $\mathbb{G}_{n}((f(U))=\sqrt{n}[\mathbb{P}_{n}(f(U))-P(f(U))]$,
$\mathbb{P}_{n}f(U)=n^{-1}\sum_{i=1}^{n}f(u_{i})$ is the empirical measure,
$Pf(U)=E(f(U))$ is the expectation taken over $U$. By empirical process
methods,
$\mathbb{P}_{n}[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)}]\rightarrow
E[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)}]=1$, in probability. Using
Slutsky’s theorem, and the empirical average converges to its limiting value.

Next,

$B_{2,m}=\frac{E[\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)}]\mathbb{G}_{n}[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)}]}{E[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)}]\mathbb{P}_{n}[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)}]}=V_{0}(\hat{d}_{j,n})\mathbb{G}_{n}(\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)})+o_{P}(1),$

by the following steps based on empirical process methods and Slutsky’s
theorem:

$\displaystyle\frac{E[\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)}]}{E[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)}]}=\frac{E[\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)}]}{E[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)}]}+o_{P}(1)=V_{0}(\hat{d}_{j,n})+o_{P}(1),$
$\displaystyle\mathbb{P}_{n}[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)}]=1+o_{P}(1),$
$\displaystyle\mathbb{G}_{n}[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)}]=\mathbb{G}_{n}[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)}]+o_{P}(1).$

Next,

$B_{3,m}=\sqrt{m}\left(\frac{E[\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X))}]-E[\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X))}]}{E[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X))}]}\right)$
$=\frac{E[Y1\\{A=\hat{d}_{j,n}(X)\\}]E[\sqrt{m}(\frac{1}{\hat{\pi}_{n}(A|X)}-\frac{1}{\pi_{0}(A|X)})]}{E[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)}]}$

$=-E[\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}\sqrt{m}(\hat{\pi}_{n}(A|X)-\pi_{0}(A|X))}{\pi^{2}_{0}(A|X)}]+o_{P}(1).$

Since
$\sqrt{n}(\pi(a|x,\hat{\theta}_{n})-\pi(a|x,\theta_{0}))=\sqrt{n}(\hat{\theta}_{n}-\theta_{0})^{T}\phi_{0}(a,x)+o_{p}(1)$,
we now have that

$B_{3,m}=-\sqrt{m/n}*\sqrt{n}(\hat{\theta}_{n}-\theta_{0})^{T}E[\phi_{0}(A,X)\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}^{2}(A|X)}]+o_{P}(1).$

Finally,

$B_{4,m}=-\sqrt{m}\frac{E(\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)})E[1\\{A=\hat{d}_{j,n}(X)\\}(\frac{1}{\hat{\pi}_{n}(A|X)}-\frac{1}{\pi_{0}(A|X)})]}{E(\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)})E(\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\hat{\pi}_{n}(A|X)})}$

$=-V_{0}(\hat{d}_{j,n})E[1\\{A=\hat{d}_{j,n}(X)\\}\sqrt{m}(\frac{1}{\hat{\pi}_{n}(A|X)}-\frac{1}{\pi_{0}(A|X)})]+o_{P}(1)$

$=V_{0}(\hat{d}_{j,n})E[\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}^{2}(A|X)}\sqrt{m}(\hat{\pi}_{n}(A|X)-\pi_{0}(A|X))]+o_{P}(1)$
$=\sqrt{m/n}\sqrt{n}(\hat{\theta}_{n}-\theta_{0})^{T}E[\phi_{0}(A,X)\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}^{2}(A|X)}]V_{0}(\hat{d}_{j,n})+o_{P}(1).$

Thus:

$\displaystyle\sqrt{m}(\hat{V}_{m}(\hat{d}_{j,n})-V_{0}(\hat{d}_{j,n}))$
$\displaystyle=\mathbb{G}_{m}(\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)})$
$\displaystyle-
V_{0}(\hat{d}_{j,n})\mathbb{G}_{m}(\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)})$
$\displaystyle-\sqrt{m/n}\sqrt{n}(\hat{\theta}_{n}-\theta_{0})^{T}E[\phi_{0}(A|X)\frac{Y1\\{A=\hat{d}_{j,n}(X)\\}}{\pi^{2}_{0}(A|X)}]$
$\displaystyle+\sqrt{m/n}\sqrt{n}(\hat{\theta}_{n}-\theta_{0})^{T}E[\phi_{0}(A|X)\frac{1\\{A=\hat{d}_{j,n}(X)\\}}{\pi^{2}_{0}(A|X)}]V_{0}(\hat{d}_{j,n})+o_{P}(1)$
$\displaystyle=\mathbb{G}_{m}(\frac{(Y-V_{0}(\hat{d}_{j,n}))1\\{A=\hat{d}_{j,n}(X)\\}}{\pi_{0}(A|X)})$
$\displaystyle-\sqrt{m/n}\sqrt{n}(\hat{\theta}_{n}-\theta_{0})^{T}E(\phi_{0}(A,X)\frac{(Y-V_{0}(\hat{d}_{j,n}))1\\{A=\hat{d}_{j,n}(x)\\}}{\pi^{2}_{0}(A|X)})+o_{P}(1),$

and the desired results follow. ∎

Proof of Proposition 2

###### Proof.

$\sqrt{m}(\hat{V}_{m}(\hat{d}_{1,n})-\hat{V}_{m}(\hat{d}_{2,n})-V_{0}(\hat{d}_{1,n})+V_{0}(\hat{d}_{2,n}))$

$=\mathbb{G}_{m}[\frac{(Y-V_{0}(\hat{d}_{1,n}))1\\{A=\hat{d}_{1,n}(x)\\}-(Y-V_{0}(\hat{d}_{2,n}))1\\{A=\hat{d}_{2,n}(x)\\}}{\pi_{0}(A|X)}]$

$-\sqrt{m/n}\sqrt{n}(\hat{\theta}_{n}-\theta_{0})^{T}E(\phi_{0}(A,X)\frac{(Y-V_{0}(\hat{d}_{1,n}))1\\{A=\hat{d}_{1,n}(x)\\}-(Y-V_{0}(\hat{d}_{2,n}))1\\{A=\hat{d}_{2,n}(x)\\}}{\pi_{0}^{2}(A|X)})+o_{P}(1)$

$\rightarrow N(0,T_{0}^{2}),$

in distribution, and standard arguments can now be used to show that
$\hat{T}_{m}^{2}\rightarrow T_{0}^{2}$ as $m\rightarrow\infty$, where

$\hat{T}_{m}^{2}=m^{-1}\sum_{i=1}^{m}(\hat{U}_{i,1}-\hat{U}_{i,2}-\bar{U}_{m,1}+\bar{U}_{m,2})^{2}+m(\hat{W}_{m,1}-\hat{W}_{m,2})^{T}\hat{\Sigma}_{n}(\hat{W}_{m,1}-\hat{W}_{m,2}).$

∎

##  Appendix B Additional Figures

Figure 5: The value function results for both training and test sets for five
different treatment regimes: observed treatment, treatment $A=0$ for all
individuals, treatment $A=1$ for all individuals, Q-learning optimal ITR, and
D-learning optimal ITR. For comparison purposes, we include the average
outcome as a benchmark. The results on original data (a) and three simulated
data sets with three levels of treatment effect modification are shown (b)
$\delta=1$, (c) $\delta=2$, (d) $\delta=10$.  Figure 6: The differences of
value functions between the pairs of the following treatments for both
training and test sets: treatment $A=0$ for all individuals, treatment $A=1$
for all individuals, Q-learning optimal ITR, and D-learning optimal ITR. The
results on original data (a) and three simulated data with three levels of
treatment effect modification are shown (b) $\delta=1$, (c) $\delta=2$, (d)
$\delta=10$.